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Abstract 



We present a general semiclassical theory of the orbital magnetic response 
of noninteracting electrons confined in two-dimensional potentials. We calcu- 
late the magnetic susceptibility of singly-connected and the persistent currents 
of multiply-connected geometries. We concentrate on the geometric effects 
by studying confinement by perfect (disorder free) potentials stressing the 
importance of the underlying classical dynamics. We demonstrate that in a 
constrained geometry the standard Landau diamagnetic response is always 
present, but is dominated by finite-size corrections of a quasi-random sign 
which may be orders of magnitude larger. These corrections are very sen- 
sitive to the nature of the classical dynamics. Systems which are integrable 
at zero magnetic field exhibit larger magnetic response than those which are 
chaotic. This difference arises from the large oscillations of the density of 
states in integrable systems due to the existence of families of periodic orbits. 
The connection between quantum and classical behavior naturally arises from 
the use of semiclassical expansions. This key tool becomes particularly sim- 
ple and insightful at finite temperature, where only short classical trajectories 
need to be kept in the expansion. In addition to the general theory for inte- 
grable systems, we analyze in detail a few typical examples of experimental 
relevance: circles, rings and square billiards. In the latter, extensive nu- 
merical calculations are used as a check for the success of the semiclassical 
analysis. We study the weak-field regime where classical trajectories remain 
essentially unaffected, the intermediate field regime where we identify new os- 
cillations characteristic for ballistic mesoscopic structures, and the high-field 
regime where the typical de Haas-van Alphen oscillations exhibit finite-size 
corrections. We address the comparison with experimental data obtained 
in high-mobility semiconductor microstriictTircs discussing the differences be- 
tween individual and ensemble measurements, and the applicability of the 
present model. 

03.65.Sq,05.45.-hb,05.30.Fk,73.20.Dx 



2 



Typeset using REVI^ 



Contents 



vr 



introduction 



A Historical perspective 



B Susceptibility of unconstrained and constrained electron systems 
C Overview of this work 



II Thermodynamic formalism 



A Grand canonical properties] 
B Canonical ensemble| . . . . 



Ill Landau susceptibilityl 



IV Systems integrable at arbitrary fields 



A Oscillating density of states for weak held 



B Circular billiards] 



C Ringi 



V Simple Regular Cjleometries: the Square 



A Oscillating density of states for weak held 



B 'I'he susceptibility: individual samples vs. ensemble average; 



C Contribution of longer orbit; 



D Numerical calculation; 



E Comparison between numerical and semiclassical results 



F Comparison with the experiment 



Cj^eneric integrable and chaotic systems 

A Generic integrable system^ 



1 Perturbation theory for magnetic helds 





2 Magnetic susceptibility for a 


generic integrable system 


B 




Generic chaotic systems 




C 


integrable vs. chaotic geometries 








I Paramagnetic character of the average susceptibility 





2 'typical magnitude of the magnetic susceptibility 



Vll 


JNon— perturbative fields: bouncing— ball— and de Haas— van Alphen- 




oscillations] 






A Intermediate helds: Bouncing- ball magnetism 






B Strong held regime 





VIII Conclusion! 



APPIINDIXES 



A Convolution of 


a rapidly oscillating function with the derivative of the 


Fermi function 





B Semiclassical expansion of the mean density of states 



73 



Calculation of for a ring billiard 



75 



D Calculation of the determinant Z^m zero field for a generic integrable 



system] 



76 





Diagonal part of the Green function for a 


free electron 


in a constant 


magnetic field] 




a, Computation of the prefactor of a Green 


function near a 


focal point 




|d Application to the cyclotron motion . . 







77 

77 



4 



I. INTRODUCTION 



A. Historical perspective 

The study of orbital magnetism in an electron gas goes back to the 1930 's with the 
pioneering work of Landau demonstrating the existence of a small diamagnetic response 
at weak fields H and low temperatures T (such that kgT exceeds the typical spacing fioj, 
uj = eH/mc). Three features of this original proposal contributed to the slowness of its 
general acceptance. First, it deals with a purely quantum result that can be expressed 
as a thermodynamic relationship without an explicit h dependence. In contrast to that 
the Bohr-van Leeuwen theorem ||^ establishes the absence of magnetism for a system of 
classical particles. For finite systems the boundary currents are shown to exactly cancel the 
diamagnetic contribution from cyclotron orbits of the interior. (This result remains valid 
even if we consider Fermi or Bose statistics 0.) Secondly, boundary effects (so crucial in 
obtaining the correct classical behavior) did not enter into Landau's original derivation. 
Twenty years later Sondheimer and Wilson ^ presented a more rigorous formulation for 
the magnetism of unconstrained electrons at weak and strong fields without using explicit 
knowledge of the energy levels, thus avoiding complicated arguments involving boundary 
electrons. (Here we present a semiclassical derivation of Landau diamagnetism independent 
of the energy level structure and valid for constrained geometries at arbitrary magnetic 
fields.) Finally, Landau diamagnetism for standard metals yields a small effect (one third of 
the Pauli spin paramagnetism) making its experimental observation rather difficult. Peierls 
showed shortly after Landau's work that the diamagnetic susceptibility persists when 
electrons are placed in a periodic potential and its value is obtained by simply using the 
effective mass instead of the free electron mass. But even if the effective mass is smaller than 
the bare mass, and the diamagnetic orbital response dominates over the spin paramagnetic 
susceptibility (as typically happens in doped semiconductors), the detailed comparison with 
the experimental data on metals was still difficult This follows from the complicated 
electronic structure and the fact that taking into account electron-electron interactions in 
the same way as a periodic potential by renormalising the effective mass, is a too crude 
approximation. 

While the restriction of the electron gas to a two-dimensional plane (still in the ther- 
modynamic limit) does not pose any new conceptual or calculational difficulty the 
effect of confining the electron system to a finite volume introduces a new energy scale in 
the problem (the typical level spacing A) and leads to a modification of the Landau suscep- 
tibility. The latter point has therefore been the object of several confiicting studie^. The 
investigation of finite-size corrections was motivated by experiments on small metal clusters 
and dealt with various theoretical models: thin plates 0, thin cylinders |]10[, confinement 



by quadratic potentials [0,0 , circular [|T3l and rectangular boxes [|14| , |T5| . Finite-size effects 



and corrections to bulk magnetism obviously depend on the relation between the typical size 
a of the system and other relevant length scales |jT6|: The thermal length Lj- = hvp/S/ir, {vp 



^For the historical account of this tortuous and often contradictory chain of findings see Ref. 
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is the Fermi velocity and j3 = 1//cbT), the elastic mean free path / (with respect to impurity 
scattering), and the phase-coherence length L$ (taking into account inelastic processes like 
electron-phonon scattering). Most of the above mentioned studies neglect other scattering 
mechanisms than that by the boundaries of the device, and deal with the macroscopic (or 
high temperature) case of Lt <^ a. The first assumption severely limits the possible com- 
parison with the experimental data of real metal clusters, while in the regime where the 
second assumption is valid the magnetic response is dominated by its smooth component 



for which, as generally shown by Robnik and Antoine |]T8|, only small corrections to the 
the diamagnetic bulk susceptibility are found. 

Opposite to the macroscopic limit, there had been studies in the extreme quantum limit 
il9| , where the temperature is low enough to enable the resolution of individual levels {kgT < 
A). In this regime the magnetic susceptibility is dominated by erratic fluctuations that, 
as we will see later, hinder its unequivocal determination. The purpose of the present 
work is the study of size corrections in the mesoscopic regime intermediate between 
the two previous limits (that is, for temperatures verifying L^/a > 1 > (3 A) and where 
inelastic processes do not inhibit quantum interference effects (L$ > a). Nowadays this is 
an experimentally accessible regime receiving considerable attention due to the richness of 
its physical properties. When L$ > a > / we have the mesoscopic dijfusive regime where the 
electron motion is dominated by impurity scattering, while for L$ > / > a we enter into the 
ballistic regime where electrons are mainly scattered off the walls of a confining potential. A 
central conclusion of our work is that finite-size corrections to the magnetic susceptibility 
in the ballistic regime can be orders of magnitude larger than the bulk values. 

One of the reasons for the sustained interest of the last few years in the mesoscopic 
ballistic regime is the possibility of studying the relation of the underlying classical dynam- 
ics and the quantum properties. This issue is precisely the subject of the field known as 



"quantum chaos" pi| , |22[| . Since the number of electrons in a mesoscopic system is always 
large, particles at the Fermi energy have a De Broglie wave length much smaller than the 
typical size a of the system {a/Xp (x kpa oc N^^"^, kp is the Fermi wave vector, d the number 
of degrees of freedom), and are therefore well in the semiclassical regime. High-mobility 
mesoscopic semiconductor samples provide an appropriate experimental testing ground in 
this context and have been recently examined with respect to the role of chaos in transport 
phenomena (for a review, see p3| , |2^ ). The present work extends the connection between 
mesoscopic systems and quantum chaos to thermodynamic properties, and analyzes recent 
experiments P5| , PB| measuring the magnetic response of ballistic microstructures. One main 
concern of our work is to show that mesoscopic finite-size effects depend crucially on the 
classical dynamics of the ballistic billiard, i.e. whether it is integrable or chaotic, and that 
the magnetic response provides an experimentally accessible criterion in order to distin- 
guish between integrable and chaotic devices (much more neatly than through the subtle 
differences found in the transport problem p3| , p^ ). 

The importance of geometrical effects for the finite-size corrections in the above defined 
macroscopic limit had already been noticed in terms of the sensitivity of the magnetic 
susceptibility on the structure of the confining potential [[10| -[T^ . The chosen potentials were 
obviously non generic but used due to their calculational simplicity, and therefore it was 
not possible to anticipate the order-of-magnitude effect that classical dynamics might have 
on the susceptibility outside the macroscopic limit. The problem of orbital magnetism from 
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a quantum chaos point of view was first addressed by Nakamura and Thomas |^ in their 



numerical study of the differences in the magnetic response of circular and elliptic billiards 
at zero temperature. The circular billiard is integrable at arbitrary field, while the ellipse 
develops chaotic behavior at finite fields. They found a reduction compared to the bulk 
susceptibility and strong fluctuations (with varying magnetic field) , and observed that both 
effects were stronger for the elliptic billiard. As already mentioned, the difficulty of these 
studies in the extreme quantum limit (at zero temperature) consists in the existence of strong 
fluctuations arising from exact or quasi-crossings of energy levels (depending parametrically 
on the magnetic field) where the susceptibility diverges. Similar features were obtained for 
other integrable systems in the quantum limit like the rectangular box the Corbino disk 
and the cylinder [^. However, this unphysical behavior is regularized by finite temperature 
that approximately adjusts the populations of both levels to each other at a crossing (or 
anti-crossing). 

Parallel to the studies of the orbital response in finite size singly-connected systems, 
there have been important developments in the understanding of persistent currents (i.e. 
the orbital magnetism in multiply-connected geometries)[|. These latter studies started 
usually from very general considerations without making the connection with the Landau 
diamagnetism. The pioneering work of Biittiker, Imry and Landauer demonstrating that 



in the presence of magnetic flux the ground state of a one dimensional ring has a current flow 
generated a large theoretical activity, mainly directed towards generalizations of quasi-one 



dimensional and diffusive rings p0| , |32| . The first experimental evidence of persistent currents 



in an ensemble of mesoscopic copper rings was given by the 1990 measurement of Levy et 



al 1^31 . The use of an ensemble was motivated by experimental reasons and brought up 
important issues about the differences between the canonical and grand canonical ensembles 
in the mesoscopic regime [pi-|5B| that we will review in the present work. Later experiments 



achieved the measurement of persistent currents in single disordered and ballistic ||2^ 



rings. In Sec. |V] we analyze in detail the last experiment making the connection with the 
orbital magnetism of the other sections. The connection between classical mechanics and 



persistent currents has already been explored in Refs. [p8H40[ . 

Small metallic samples at sufficiently low temperatures operate in the diffusive meso- 
scopic regime, where the classical electron motion is a random walk through the impurity 
potential. This was the regime of the original experiment on persistent currents [^] and 
therefore received considerable theoretical attention. The effect of disorder on persistent 



currents has been evaluated by diagrammatic perturbation theory |]5B|. A weak disorder 



potential does not alter the bulk Landau diamagnetism and gives within perturbation 



theory enhancement factors proportional to kpl for finite samples Highly pure semi- 
conductor heterojunctions combined with lithographic techniques allow the realization of 
samples small enough to be in the mesoscopic ballistic regime where I > a. This is the 



case of the orbital magnetism and persistent current measurements of Refs. and ^6 
In the ballistic regime electrons move almost straight between collisions with the walls of 
the confining potential. The small drift between collisions is due to the unavoidable disor- 



^For historical accounts on persistent currents see Refs. |3C,31|. 
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der potential existing in real structures. Neglecting completely the effect of disorder, and 
therefore the associated drift, leads to an ideal or dean system which describes simply an 
electron billiard. A central result of our work is that the application of semiclassical expan- 
sions at finite temperature allows one only to consider short classical periodic trajectories. 
Therefore, the clean model provides a reasonable approach to the weak and smooth disorder 
of the ballistic regime. We will get back to this point in this work, and in a separate paper 
PB| we examine in detail the role of disorder in ballistic samples. 

The previously cited developments, as well as most of the present work, deal with finite- 
size effects in the orbital response at weak fields. At high fields the magnetic response is 
dominated by the occurrence of Landau levels (whose spacing hu is much larger than A 
or kgT) yielding the well known de Haas-van Alphen effect. In 1938 Landau derived (see 
Refs. and a complete analytical expression for the susceptibility of a degenerate free 
electron gas including the weak-field diamagnetic response and the de Haas-van Alphen 
oscillations. Since the latter turned out to be a powerful technique to examine the elec- 
tron structure of metals its study has been at the heart of Condensed Matter Physics 
for various decades. Its measurement in two-dimensional electron gases has allowed the 



determination of the density of states at high field With the advent of Mesoscopic 

Physics the question of finite-size effects on the de Haas - van Alphen effect was naturally 
raised, and free electron gases on a disk |46,|^ and confined by a parabolic potential | i8[ 



were considered at high fields. The semiclassical theory used in the present work provides 
the finite-temperature susceptibility at arbitrary fields and allows the identification of an 



intermediate regime characteristic for ballistic samples that we discuss in Sec. [VII 



Our work aims at the convergence of various seemingly disconnected fields: Landau 
diamagnetism, persistent currents, de Haas - van Alphen effect, finite-size corrections of 
thermodynamic functions, quantum chaos, and electronic properties of weakly disordered 
systems. We will show that the semiclassical analysis naturally enters in the problem of the 
magnetic response of ballistic structures, provided a model of non-interacting electrons is 
adequate. The expression of the magnetic susceptibility and persistent currents in terms of 
classical trajectories provides a unifying approach applicable to various geometrical shapes, 
different temperatures and magnetic field strengths. 



B. Susceptibility of unconstrained and constrained electron systems 

We now present the basic formulas defining the magnetic susceptibility and then com- 
pare the unconstrained magnetic response with the susceptibility obtained by confining the 
electron gas to a finite region to illustrate the subject of our studies in this paper. Let us 
consider a noninteracting electron gas confined in a volume (area in two dimensions) A at 
temperature T under a magnetic field H. The magnetic moment of the system in statistical 
equilibrium is given by the thermodynamic relation 



-^ = -^77 (1-1) 




where Q{T, n, H) is the thermodynamic potential, and /i the chemical potential of the elec- 
tron gas. The differential magnetic susceptibility is defined by 
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dip 



(1.2) 



The notation with the superscript GC is used in order to emphasize the fact that we are 
working in the grand canonical ensemble. The choice of the ensemble in the macroscopic 
limit of N and A — oo is a matter of convenience. As it is well known by now P3-PB| the 



equivalence between the ensembles may break down in the mesoscopic regime that interests 
us, and this point will be thoroughly discussed in the remaining of the paper. However, for 
the purpose of this didactical introduction we will work in the grand canonical ensemble 
studying the magnetic response of electron systems with fixed chemical potentials. The 
calculation advantages of the GC ensemble arise from the simple form of the thermodynamic 
potential 

n{T,fi,H) = -^J dE d{E) ln(l + exp[/3(/i-E)]) , (1.3) 

in terms of the single-particle density of states 

d{E) = g,Y.^{E - E,) . (1.4) 

A 

The factor gs = 2 takes into account spin degeneracy. Ex are the eigenenergies of the system. 
The magnetic susceptibility is directly extracted from the knowledge of the density of states. 
The case of a free electron gas is particularly simple since the electron eigenstates are Landau 
states with energies 

Ek = hw{k + 1/2) fc = 0,1,2,... (1.5) 

and degeneracies gs$/$o- The cyclotron frequency is w = eH/mc, $ = HA is the flux 
through an area A, and $0 = hc/e is the elemental flux quantum. Throughout this work 
we will neglect the Zeeman splitting term due to the electron spin. It can however be 
incorporated easily when spin-orbit coupling is negligible 

Landau's derivation of the magnetic susceptibility of a free electron system arising from 
the quantization condition ( [1.5|) can be found for the three-dimensional case in standard 
textbooks The two-dimensional case [|^,|^ follows upon the same lines. In the following 
we present a sketch of the latter which will be useful towards a semiclassical understanding 
of the problem. {H is now the component of the field perpendicular to the plane of the 
electrons.) 

By the use of the Poisson summation formula the density of states related to the quan- 



tization condition (1.5) can be written as 



mA mA^, /27mE\ , , 



This decomposition is usually interpreted as coming from the Weyl term (given by the volume 
of the energy manifold in phase space) and the contribution of cyclotron orbits (second term, 
strongly energy dependent). We stress though that in the bottom of the spectra, from which 
the Landau diamagnetic component originates, this distinction is essentially meaningless. 
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In the case of a degenerate electron gas with a weak field such that hw ^ ksT ^ /i the 
energy integral ( p. .31) is easily performed resulting in 

n(fi) ~ = -gs^ ^ + gs^—^ ^ , (1.7) 

where is the smooth part (in energy) of the thermodynamic potential. (Note that the 
second term of Eq. ( |1.7| ) comes nevertheless from the integral of the rapidly oscillating 
term of the density of states.) Thus, we obtain the two-dimensional diamagnetic Landau 
susceptibility 

For high magnetic fields, ksT <^ hw, the energy integrals are slightly more complicated 
than before since the rapidly oscillating component of Q is not negligible any longer. This 
latter can be computed (see Appendix |^ for the treatment of similar cases) as 

n---.{^)U-n^(^J^^^)Mn). (1.9, 

where -Rt(^) is a temperature dependent damping factor 

27r'^nkBT/hw 

^^^^^ = sinh(2vr%A:^r/;^.^) " ^^'^^^ 

With = + we have the Landau and de Haas-van Alphen contributions to the 
magnetic susceptibility 

-i-24(£)=E(-irco.(^)..,«). (i.n, 



Xl 



n=l 



The second term exhibits the characteristic oscillations with period 1/H and is exponentially 
damped with temperature (and the summation index n)^ 

While going from the bulk two-dimensional case (macroscopic regime) to the constrained 
case (ballistic mesoscopic) two important changes take place: i) the confining energy appears 
as a relevant scale and Eq. ( |1.5| ) no longer provides the quantization condition; ii) since we 
are not in the thermodynamic limit of and A ^ oo, the constraint of a constant number 
of electrons in [isolated] microstructures is no longer equivalent to having a fixed chemical 
potential. These two effects will be thoroughly discussed in the paper. For didactical 
purposes we restrict ourselves in this introductory section to only the changes (i) due to the 
confinement, and we anticipate some of the results that will be later discussed in detail. 



^For high fields we cannot in principle separate the orbital and spin effects. The de Haas-van 
Alphen oscillations are given only by the orbital component, that is the only one that interests us 
for our model of spinless electrons. 
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We imagine a mesoscopic square of size a connected to an electron reservoir with chemical 
potential fi. Direct numerical diagonalization in the presence of a magnetic field (Fig. |l]a) 
allows us to obtain x^'^ (solid line in Fig. |I]b). In the high field region (2rc < a, we note 
fc = Vp/uj the cyclotron radius) the characteristic de Haas - van Alphen oscillations are ob- 
tained, although not with the amplitude expected from calculations in the bulk (Eq. ( |1 . 11| ) ) ■ 
For lower fields the discrepancy between our numerical results and the bulk Landau diamag- 
netism is quite striking. Thus, confining deeply alters the orbital response of an electron 
gas. Without entering into details at this point we remark the fact that the whole curve is 
quite well reproduced by a finite-temperature semiclassical theory (dashed line) that takes 
into account only one type of trajectory (see insets) in each of the three regimes: a) the 
interference-like regime, dominated by the shortest trajectories with the largest enclosed 
area for squares at zero magnetic field; b) the transition regime dominated by the bend- 
ing of bouncing-ball trajectories between parallel sides of the square; c) the de Haas - van 
Alphen regime dominated by cyclotron orbits. It is remarkable how an exceedingly com- 
phcated spectrum as that of Fig. |l|a can be understood within such a simple semiclassical 
picture once finite temperature acts as a filter selecting only few types of trajectories. 



C. Overview of this work 

The purpose of this paper is to provide an [essentially comprehensive] theory of the orbital 
magnetic properties of non-interacting spinless electrons in the mesoscopic ballistic regime. 
We restrict ourselves to the clean limit, where the different behavior of the magnetic response 
arises as a geometrical effect (shape of the microstructure). We will make extensive use of 
semiclassical techniques since they appear to be perfectly suited for these problems. For 
the smooth components (such as in Eq. (|1.7|)) we will use the general techniques developed 
by Wigner to obtain higher h corrections to the Weyl term which are field dependent. For 
the oscillating components, we will rely on the so called semiclassical trace formulas, which 
provide simple and intuitive expressions for the density of states as a sum over Fourier-like 
components associated to closed classical orbits. 

In this respect it will be seen that the nature of the classical dynamics, i.e. integrable 
versus chaotic (and more precisely existence versus absence of continuous families of periodic 
orbits), plays a major role. Although we will present a complete formalism for both cases, our 
main emphasis, and in particular all the examples treated explicitly, will concern integrable 
geometries. The reason for this choice is twofold. First, as we will make clear in the sequel, 
one expects a much larger magnetic response for integrable systems than for chaotic ones, 
yielding a more striking effect easier to observe. The second point is that, contrarily to 
what might seem natural a priori, integrable geometries present a few conceptual difficulties 
in their treatment which are not present for chaotic systems. Indeed integrable systems 
lack of structural stability, which means that under any small perturbation (such as the 
one provided by the presence of a magnetic field) they generically do not remain integrable. 
Chaotic systems on the contrary remain chaotic under a small perturbation. Therefore, 
as shown in Ref. [Q, the Gutzwiller trace formula pl| , |5l| , valid for chaotic systems, can 



be used at finite fields without further complications. For integrable geometries however. 



the Berry- Tabor [p2| , ^ or Balian-Bloch trace formulae valid for integrable systems 
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usually do not apply in the presence of a perturbing magnetic field. It will therefore be 
necessary, following Ozorio de Almeida to consider the more complicated case of 

nearly integrable systems, which we will do in detail here. 

To perform this program, the present work is organized as follows. In the next section 
we present the thermodynamic formalism appropriate for working in the canonical and 
grand canonical ensembles, stressing its semiclassical interpretation and incorporating the 
changes due to the constancy of the number of electrons in the experimentally relevant 
microstructures. In Sec. |T| we consider the smooth magnetic response and show that the 
Landau diamagnetism is present in any confined geometry at arbitrary fields. In Sec. |TV| we 
address the magnetic response (susceptibility and persistent currents) in the simplest possible 
geometries: circles and rings billiards that are integrable with and without magnetic field. 
In Sec. |V] we present the calculation of the magnetic susceptibility for the experimentally 
relevant case of the square billiard whose integrability at zero field is broken by the 
effect of an applied magnetic field. An initial study along these lines was presented in 
Refs. and and independently proposed by von Oppen |^. This geometry and 



the corresponding experiment have also been analyzed from a completely different point 
of view by Gefen, Braun and Montambaux stressing the importance of the residual 
disorder (see also Ref. [§![]). We consider in Sec. |VI| the generic magnetic response of both 
integrable and chaotic geometries, stressing the similarities and differences in their behavior 
and calculating the line-shape of the average magnetization in generic chaotic systems. In 
Sec. [VI I| we demonstrate how the semiclassical formalism we have developed applies not only 
to the weak-field limit, but also to higher field and in particular to the high field regime of the 
de Haas - van Alphen oscillations. We treat explicitly the example of the square geometry, 
including an intermediate field regime dominated by bouncing-ball orbits as depicted in 
Fig. We discuss our conclusions and their experimental relevance in Sec. |V111| . The 
modifications of our results due to the effect of a weak disordered potential are discussed in 
a separate publication 

To keep the focus on the physical concepts developed in the text, a few technical deriva- 
tions have been relegated to some appendices. Appendix ^ presents the generic case of 
the convolution of a rapidly oscillating function with the derivative of the Fermi function. 
Appendix ^ gives the calculation of the first field-dependent term of the heat Kernel in an h 
expansion. In Appendix |C| we compute the action integrals associated with the dynamics of 
circular and ring billiards needed to define the energy manifold in action space. Appendix 
presents the calculation of the prefactor of the Green function for an integrable system, 
while in Appendix ^ we show how to compute the semiclassical Green function at a focal 
point, and apply the obtained result to the particular case of cyclotron motion. 
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II. THERMODYNAMIC FORMALISM 



One main subject of the present work is the introduction of semiclassical concepts into 
the thermodynamics of mesoscopic systems. In this section we provide the basic formahsm 
allowing one to obtain the thermodynamic properties (grand potential, free energy) from 
the quasiclassically calculated single-particle density of states and hence the susceptibility. 
We begin with general definitions and relations between grand canonical and canonical 
quantities. 

For a system of electrons in a volume (area in two dimensions) A connected to a reservoir 
of particles with chemical potential /x (grand canonical ensemble) the magnetic susceptibility 
is obtained, as given by Eq. (|1.2|) , as 



fl{T, fi, H) is the thermodynamic potential, which can be expressed for non-interacting 
electrons in terms of the single-particle density of states through Eq. ( [1.3|) . 

For actual microstructures, the number N of particles inside the device might be large but 
is fixed in contrast to the chemical potential /i. As discussed in the introduction, it will be 
necessary in some cases, namely when considering the average susceptibility of an ensemble 
of microstructures, to take explicitly into account the conservation of N, and to work within 
the canonical ensemble. For such systems with a fixed number N of particles, the relevant 
thermodynamic function is not the grand potential Q, but its Legendre transform, the free 
energy^ 



Except for the calculation of the Landau contribution performed in the following section 
all the computations of the magnetic response of the microstructures to be considered will 
involve two clearly separated parts. In the first one the (oscillating part of the) density of 
states will be calculated semiclassically. Depending on the underlying classical dynamics 
(integrable versus chaotic, with or without breaking of the invariant tori, with or without 



^In standard thermodynamics, Eq. ( |2.1[ ) just represents the definition of the grand potential. 
It should be borne in mind however that from a statistical physics point of view this is not an 
exact relation, but the result of a stationary-phase evaluation of the average over the occupation 
number, valid only when ksT is larger than the typical level spacing. Therefore, we are entitled 
to use this relation in the mesoscopic regime that interests us, but not in the microscopic regime, 
where features on the scale of a mean spacing become relevant. 






In particular, the magnetic susceptibility of a system of N electrons is 




(2.2) 
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focal points, etc.), the results as well as their derivation will vary noticeably. In the second 
stage the integrals over energy yielding the desired thermodynamic properties have to be 
performed in a leading order in h approximation. To avoid tedious repetitions, we shall 
consider here in some detail this part of the calculation of the thermodynamic properties, 
and refer without many additional comments to the results obtained in this section whenever 
needed. We begin with the grand canonical quantities which exhibit the simplest expressions 
in terms of the density of states. In a second subsection we shall consider the canonical 
ensemble following closely the approaches presented in Refs. . 



A. Grand canonical properties 

We begin with the standard definition, Eq. (|1.4|) of the density of states 



d{E) = g^Y^^^E - Ey) , 

A 

(gs = 2 is the spin degeneracy. Ex the eigenenergies) and its successive energy integrals. 
They are the energy staircase 

niE) = dE' d(E') , (2.3) 
Jo 

and the grand potential at zero temperature 

uj{E) = - ['^ dE' n{E') . (2.4) 
Jo 

These are purely quantum mechanical quantities, depending only on the eigenstates Ex of 
the system. At finite temperature the corresponding quantities are obtained by convolution 
with the derivative f'{E — n) of the Fermi distribution function 



We then have 



fOO 

D{fi) = - dE d{E) f'iE-fi) , (2.6a) 
Jo 

POO 

N{fi) = - dE n{E) f'{E-^i) , (2.6b) 
Jo 

POO 

Q{fi) = - / dE uo{E) f'iE-fi) . (2.6c) 



Integration by parts leads to the standard definition (|1.3|) of the grand potential and the 
mean number of particles in the GCE with a chemical potential /x, i.e. 

POO 

N{fi) = / dE d{E) f{E-fi) . (2.7) 
Jo 



That means that the thermodynamic properties (|2.6b|) - (|2.6c|) are obtained by performing 



the energy integrations (|2.3|)-(|27^) with the Fermi function as a weighting factor. 
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In the following the separation of the above quantum mechanical and thermodynamic 
expressions into smooth (noted with a " " ") and oscillating (noted with the superscript 
" ") parts is going to play a major role. It has its origin in the well-known decomposition 
of the density of states as 

d{E) = d{E) + rf°^=(E) . (2.8) 

This decomposition has a rigorous meaning only in the semiclassical {E — >■ oo) regime 
for which the scales of variation of d and d°^'^ decouple. To leading order in h, the mean 
component d{E) is the Weyl term reflecting the volume of accessible classical phase space 
at energy E (zero-length trajectories), while d°'^'^{E) is given as a sum over periodic orbits 
(Gutzwiller and Berry- Tabor trace formulas) Ell]. Generically, it will be expressed as a sum 



^osc^^) = ^ dtiE) ; dtiE) = AiE) sin iSt{E)/h + u^) . (2.9) 
t 

running over periodic orbits labeled by t where St is the action integral along the orbit t, 
At{E) is a slowly varying prefactor and Vt a constant phase.0 

Using the expression ( |2.9| ) for d"'^'^ in Eqs. ( |2.^ ) and ( |2.4| ), n"'^'^ and uj°^'^ are obtained to 
leading order in % as 

7f^^{E) = dE' d^'^^iE') ; cu^'^iE) = - dE' n°''{E') . (2.10) 

The lower bounds are not specified because the constants of integration are determined by 
the constraint that n'^^'^ and uj°^'^ must have zero mean values. (It should be borne in mind 
that semiclassical expressions like ( |2.8D , and those that will follow, are not applicable at the 
bottom of the spectrum.) 

In a leading h calculation the integration over energy in Eq. ( |2.1CI| ) has to be applied only 
to the rapidly oscillating part of each periodic orbit contribution dt. Noting moreover that 
if St{E) is the action along a periodic orbit, then Tt{E) = dSt/dE is the period of the orbit, 
one has in a leading h approximation 

f^A{E') sm{St{E')/h + Ut)dE' = ^A{E) cos {St{E)/h + u^) (2.11) 

as can be checked by differentiating both sides of Eq. (|2.11|) . In order to emphasis that 
the integration over energy merely yields a multiplication by {—h/r), we use the notation 
{i^ ■ dt) to assign the contribution dt of a periodic orbit after shift of the phase by n/2, i.e. 
(i^ ■ [B sm{S/h)]) = B cos{S/h). We get 

n'^''iE) = Y.ntiE) ; utiE) = ^ (t^ ■ dtiE)) , (2.12) 




^osc(^) = ^cUi(E) ; UtiE) = -— dtiE) . (2.13) 



^When considering systems whose integrability is broken by a perturbing magnetic field, we shall 
stress the necessity to consider families of recurrent - but not periodic - orbits of the perturbed 
system. This will, however, not affect the discussion which follows. 
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The thermodynamic functions D"^'^{fi), N°^^{fi) and are then obtained by apph- 

cation of Eqs. (|2.6|) in which the full functions are replaced by their oscillating component. 
The resulting integrals involve the convolution of functions {d"'^'^{E), n"'^'^{E) or 00°^'^ (E)) 
oscillating [locally around /i] with a frequency t{^)/{2ttK), with the derivative of the Fermi 
factor f'{E — fi) being smooth on the scale of = k^T. One can therefore already antici- 
pate that this convolution yields an exponential damping of the periodic orbit contribution 
whenever r(/x) 3> h/S. As shown in appendix ^ the temperature smoothing gives rise to an 
additional factor for each periodic orbit contribution, 

^T{Tt) = . , . , ^ , = — , (2.14) 

smh(n/rc) TT 

in a leading h and (3~^ approximation (without any assumption concerning the order the lim- 
its are taken). In this way we obtain relations between the following useful thermodynamic 
functions and the semiclassical density of states: 

^°''(/i) = EA(/x) ; Dt{fi) = RTin)dt{fi) , (2.15a) 
t 

iV°^^(/x) = ; ^*(/^) = Mn) ■ dti^M)) , (2.15b) 

^°''(/^) = E^*(/^) ; n^iij,) = Rrin) (^^^ dtifx) . (2.15c) 

At very low temperature, Rt ~ 1 — [(r(7r)/(/i/5)]^/6 which, for billiard-like systems where 
Tt = Lt/vp (with Lt being the length of the orbit and Vp the Fermi velocity), simply gives 
the standard Sommerfeld-expansion R^ ~ 1 — [(Lt7r)/(^/3uj^)]^/6. For long trajectories or 
high temperature it yields an exponential suppression and therefore the only trajectories 
contributing significantly to the thermodynamic functions are those with < Tc- Thus, 
temperature smoothing has a noticeable effect on the oscillating quantities since it effec- 
tively suppresses the higher harmonics, which are associated with long classical orbits in a 
semiclassical treatment. On the contrary, for a degenerate electron gas (/9/x ^ 1), finite tem- 
perature has no effect on the mean quantities. Temperature is then the tuning parameter 
for passing from d{E) at T = to D[E) = d{E) at large temperatures (by the progres- 
sive reduction of d"^'^). Similar considerations hold for the energy staircase and the grand 
potential. 

The oscillatory part of the semiclassical susceptibility in the grand-canonical ensemble 
is finally obtained from Eq. ( |1.2|) by replacing Q by 



B. Canonical ensemble 

Let us now consider the susceptibility in the canonical ensemble, appropriate for systems 
with a fixed number of particles. We follow Imry's derivation for persistent currents in 
ensembles of disordered rings [^]. The only important difference is that we will take averages 
over the size and the Fermi energy of ballistic structures instead of averages over impurity 
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realizations. We will stress the semiclassical interpretation that will be at the heart of our 
work, and highlight some of its subtleties. 

As mentioned in the introduction the definition Eq. ( p.2| ) of the susceptibility x is equiv- 
alent to up to 1/N {i.e. h corrections). Therefore, in the macroscopic limit of N — >■ oo 
the choice of the ensemble in which the calculations are done is unimportant. On the other 
hand, in the mesoscopic regime of small structures (with large but finite N) we have to 
consider such corrections if we want to take advantage of the computational simplicity of 
the Grand Canonical Ensemble (GCE). The difference between the two definitions is par- 
ticularly important when the GCE result is zero as it is the case for the ensemble average of 
X*^*^. The evaluation of the corrective terms can be obtained from the relationship Eq. 
between the thermodynamic functions -F(N) and and the relation N{^) = N. In the 

case of finite systems the previous implicit relation is difficult to invert. However, when N is 
large we can use the decomposition of N{fi) in a smooth part N{^) and a small component 
N°^'^{fi) that fluctuates around the secular part, and we can perturbatively treat the previous 
implicit relation. The contribution of a given orbit to d°^^ is always of lower order in h than 
d as can be checked for the various examples we are going to consider and by inspection of 
semiclassical trace formulae. However, since there are infinitely many of such contributions, 
we obtain d°^'^ and d to be of the same order when adding them up. (This must be the case 
since the quantum mechanical d{E) is a sum of 6 peaks.) Thus, we cannot use d°^^/d as a 
small expansion parameter. On the other hand, finite temperature provides an exponential 
cutoff in the length of the trajectories contributing to D"^'^ so that only a finite number of 
them must be taken into account. Therefore, D"^'^ is of lower order in h than Z), and in 
the semiclassical regime it is possible to expand the free energy F with respect to the small 



parameter D°^'^/D. The use of a temperature smoothed density of states Eq. ( 2.6a ) closely 



follows the Balian and Bloch approach [^, where, due to the exponential proliferation of 
orbits and the impossibility of exchanging the infinite time and semiclassical limits, the semi- 
classical techniques based on trace formulae are considered meaningful only when applied 
to smoothed quantities. The decomposition of D{E) is depicted in Fig. ^, where we have 
taken Z)(~ d) to be energy independent, corresponding to the two-dimensional (potential 
free) case. 

For a perturbative treatment of the mentioned implicit relation we define a mean chemical 
potential fi by the condition of accommodating N electrons to the mean number of states 

N = N{ij) = N{fL) . (2.16) 

Expanding this relation to first order in D^'^'^/D, and employing that dN/dfi = D, one has 

A;x = /.-/x^- J-y iV°-(/.). (2.17) 

The physical interpretation of A/i is very clear from Fig. § The shaded area represents the 
number of electrons in the system and it is equal to the product D x fi. 



^In the following we will only write the N dependence of F and the ^ dependence of J7, assuming 
always the T and H dependence of both functions. 
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Expanding the relationship ( |2.1|) to second order in A/x, 



F(N) ~ (/i + A/i)N + fi(^) - N{fi)Afi - D{ji) , (2.18) 

using the decomposition of and N{fi) into mean and oscillating parts and eliminating 
A/x (Eq. ( p.l7|) ) in the second order term, one obtains the expansion of the free energy to 



second order in D°^'^/D ||35| , p6 



F(N) ~ F° + AF(^) + AF(2) , (2.19) 

with 

F° = + n{fx) , (2.20a) 

AF(i) = Q°'^{fi) , (2.20b) 

AF(2) = ^ (iV°-(/.))^ . (2.20c) 

Then AF^^) and AF^^) can be expressed in terms of the oscillating part of the density of 
states by means of Eqs. (1215^1) and (|2J[5^ ). The first two terms F° + AF^^) yield the 
magnetic response calculated in the GCE with an effective chemical potential /2. The first 
"canonical correction" AF^^^ has a grand canonical form since it is expressed in terms of 
a temperature smoothed integral of the density of states (Eq. (^T^)) for a fixed chemical 
potential fl. 

It is convenient to use the expansion ( p.l9|) in the calculation of the magnetic suscep- 
tibility of a system with a fixed number of particles because the leading h contribution to 
iV(/i) has no magnetic field dependence, independent of the precise system under considera- 
tion. Therefore, at this level of approximation, keeping N constant in Eq. (|2.2|) when taking 
the derivative with respect to the magnetic field amounts to keep /i constant. Since F^'^^ is 
field independent in a leading order semiclassical expansion the weak-field susceptibility of 
a given mesoscopic sample will be dominated by AF*^^^ . However, when considering ensem- 
bles of mesoscopic devices, with slightly different sizes or electron fillings, AF^^^ (and its 
associated contribution to the susceptibility) averages to zero due to its oscillatory behavior 
independently of the order in h up to which it is calculated.[] Then we must consider the 
next order term AF^^\ 

As mentioned above, we will essentially work in the semiclassical regime (to leading order 
in h) where F° is field independent. However, in the following section we will examine the 
next order h correction to fi{fl) (and to F°), demonstrating that its field dependence gives 
rise to the standard Landau diamagnetism, independent of any confinement. 



^In the following, we shall always calculate AF^^^ in a leading order h approximation. Higher 
order corrections to AF^^^ may be of the same order as AF^^^ but will average to zero under 
ensemble averaging. 
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III. LANDAU SUSCEPTIBILITY 



In the previous section we showed that the various quantum mechanical {i.e. d{E), n{E), 
ijj{E)) and thermodynamic {i.e. D{^), N{^), properties of a mesoscopic system can 

be decomposed into smooth and fluctuating parts. In the semiclassical limit, where the 
Fermi wavelength is much smaller than the system size, each of these quantities allows an 
asymptotic expansion in powers of k. For most of the purposes it is sufficient to consider 
only leading order terms while higher order corrections must only be added if the former 
vanish for some reason. This is the case for the smooth part 0,{n) of the grand potential, 
which is the dominant term at any temperature, but is magnetic field independent to leading 
order in h. The present section will be the only part of our work where higher h corrections 
are considered. We will show that they give rise to the standard Landau susceptibility. Our 
derivation relies neither, on the quantum side, on the existence of Landau levels, nor, on the 
classical side, on boundary trajectories or the presence of circular cyclotronic orbits fitting 
into the confinement potential. This shows that the Landau susceptibility is a property of 
mesoscopic devices as well as infinite systems, being the dominant contribution at sufficiently 
high temperature^. 

We consider a (i-dimensional {d = 2, 3) system of electrons governed by the quantum 
Hamiltonian 

n=^ fp--A(q)V + y(q), (3.1) 



2m V c 

where A is the vector potential generating the magnetic field H and V^(q) is the potential 
which confines the electrons in some region of the space. This region can a priori have 
any dimension, and it can be smaller that the cyclotron radius. We will only assume in 
the following that V{q) is smooth on the scale of a Fermi wavelength, so that semiclassical 
asymptotic results can be used. In billiards the effect of hard boundaries on the susceptibility 



is negligible compared to the Landau bulk term ||17| , P^ , and therefore the results obtained 
below apply there, too. 

There exist general techniques to compute the semiclassical expansion of the mean part 
of the density of states (or of its integrated versions Eqs. ( |2.3| ), ( p.4|) ) up to arbitrary order 
in h. The most complete approach, which allows one to take into account the effect of sharp 
boundaries, can be found in the work of Seeley |Q . However, assuming the smoothness of 



V{q) allows us to follow the standard approach introduced by Wigner in 1932 |^ which is 
based on the notion of the Wigner transform of an operator. As a starting point we consider 
the Laplace transform of the level density (or heat Kernel), 



Z{\) = I dE e-^^ d{E) = gs Tr(e-^^) , (3.2) 



where gs = 2 takes into account the spin degeneracy. In appendix we apply after a 
brief description the technique to calculate the first two terms of the expansion of Z{X) with 



^Analog results have been independently obtained by S.D. Prado et al. The Wigner distri- 
bution function was previously used by R. Kubo IHS] in the study of Landau diamagnetism. 
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respect to A. They yield under the inverse transformation the first two terms of the expansion 
of d{E) in powers of h. The oscillating part d"'^'^{E) of d{E) is not included in this procedure 
since it is associated with exponentially small terms in Z{X), that is, Z{X) ~ Z{X) for 
A ~ 0. This well known property can be easily seen from the integral treated in appendix |^ 
by identifying (3 with A and using the exponential form of the distribution function in the 
classical limit of high temperatures <^ 1). 

Noting 7i(q, p) the classical Hamiltonian corresponding to Eq. ( p. 11 ), the leading order 
[Weyl] contribution to Z{X) is given by Eq. (|B6|), 



^w(A) = ^2^^)d /dqdp exp(-A7^(q,p)) , (3.3) 
and the inverse Laplace transform yields the familiar result 

d^{E) = dy,{E) = |dqdp5(E-H(q,p)) . (3.4) 

In the above integrals, the substitution 

P ^ P' = P - -A (3.5) 

c 

eliminates any field dependence. Therefore 

cjw(^) = c^w(^) = - rdE' r dE" dyj{E") , (3.6) 

Jo JO 



as well as the leading term r2w(/^) of the grand potential (obtained in the high temperature 
limit of Eq. (|2.6c| )), are field independent. This is the reason for the absence of orbital 
magnetism in classical mechanics. To observe a field dependence, one must consider the 
first correcting term of Z{X) which, as shown in appendix (Eq. ( |B11| )), is given by 

Z,iX,H) = -X'^Z^ + Z',. (3.7) 
6 

Here, fiB = {eh)/{2mc) is the Bohr magneton, and Z^ = Zi{h=o) is a field independent term 
that we will drop from now on since it does not contribute to the susceptibility. 

The integrated functions n{E) and uj{E) can be obtained from their Laplace transforms 



Z{X) Z{X) 

— ' "^^^ = 

Then the first correction to the zero-temperature grand potential is 



n(A) = ^ , wiX) = . (3.8) 



u^{E) = u,{E) = ME) . (3.9) 

6 

After convolution with the derivative of the Fermi function (Eq. p.6c| )) we obtain the first 
corrective term of the grand potential 

fii(/i) = ni(/.) = 4^Dw(/x). (3.10) 
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In the grand canonical ensemble, the above equation readily gives the leading contribu- 
tion to the susceptibility 



r" = -||5w, (3.11) 

coming from the mean part of the grand potential. In Eq. ( p.ll| ) A is the confining volume 
(area for d = 2) of the electrons. Noting that Z)w = diVw/d/x, one recognizes the familiar 
result of Landau |1[]. For systems without potential (bulk, or billiard systems), it gives in 
the degenerate case (/?/x ^ 1) in two, respectively, three dimensions 

GC _ Sse^ GC _ gse^kp 

In the non-degenerate limit the susceptibility is 

The temperature independence in the degenerate regime and the power-law decay in the non- 
degenerate limit cause the dominance of the Landau contribution at high temperatures since, 
as mentioned in the previous section (and demonstrated in appendix the contributions 
from AF^^^ and AF^^^ (Eqs. (|2.20b ) and (|2.20c )) are exponentially damped by temperature. 



The Landau diamagnetism is usually derived for free electrons or for a quadratic confining 
potential [Q,^ . We have provided here its generalization to any confining potential (including 
systems smaller than the cyclotron radius). 

For a system with fixed number N of electrons, defining a Weyl chemical potential /^w 

by 

N = iVw(/iw) (3.14) 



and following the same procedure as in Sec. II B one can write 



F(°)(N) ~Fw + ni(/iw) , (3.15) 

where both /xw and 

Fw = /iwN + r2w(/iw) (3.16) 

are field independent. Therefore, the smooth part of the free energy gives the same con- 
tribution than Eq. ( p.ll| ): We recover the Landau diamagnetic response in the canonical 



ensemble, too. 

At the end of this section we would like to comment on the case of free electrons in 
two dimensions. Since Eq. ( |1.6| ) represents an exact formula for the density of states, 
d{E) = [gsfnA] / can be interpreted as the exact mean density of states, and 
^osc(^^) = (g,m/l)/(7r;i^)E^=i(-l)"cos((27rnF)/(^w)) as the exact oscillating part. How- 
ever, uj{E) being obtained by integrating d{E) twice, has a mean value which, in addition 
to —dE'^/2, contains the term (/i^if^/6)(i yielding the Landau susceptibility. In the usual 
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derivation, this term comes from the integration of d°^'^{E), more precisely from the bound- 
ary contribution at E = (i.e., from levels too close to the ground state in order to properly 
separate the mean value from oscillating parts). One should be aware that u{E) cannot be 
defined by Eq. ( |3.6|) as soon as non leading terms are considered. For this reason some care 
was required for the definitions of the last section (see the discussion around Eqs. (|2.10| )- 

(ill))- 
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IV. SYSTEMS INTEGRABLE AT ARBITRARY FIELDS 



In the remainder of this work we will provide semiclassical approximations for the cor- 
rective free-energy terms AF(i) and AF^^) (see Eq. ( pTTOD ) and their associated magnetic 
responses for systems that react differently under the influence of an applied field. We will 
be mainly working in the weak- field regime (except in section |VI1|) , where the magnetic field 
acts as a perturbation almost without altering the classical dynamics. In this regime the 
nature of the zero-field dynamics {i.e. integrable vs. chaotic, or more precisely, the orga- 
nization of periodic orbits in phase space) becomes the dominant factor determining the 
behavior and magnitude of the magnetic susceptibility. For systems which are integrable 
at zero field the generic situation is that the magnetic field breaks the integrability (as any 
perturbation will do). It is necessary in that case to develop semiclassical methods allowing 
to deal with nearly, but not exactly, integrable systems. This question will be addressed in 
sections [V] and |V^. There exist however "non generic" systems where the classical dynamics 
remains integrable in the presence of the magnetic field. Due to their rotational symme- 
try, circles and rings (which are the geometries used in many experiments) fall into this 
category. In these cases (and similarly for the Bohm-Aharonov flux |3^) the Berry- Tabor 
semiclassical trace formula []52| , |53| provides the appropriate path to calculate semiclassically 
the oscillating part of the density of states d°^'^, including its field dependence. Thus, AF^^^ 
and AF^'^\ and their respective contributions to the susceptibihty, can be deduced. This 
is the program we perform in this section, treating specifically the example of circular and 
ring billiards. 

The magnetic susceptibility of the circular billiard can be calculated from its exact quan- 
tum mechanical solution in terms of Bessel functions 1 10 1^,30]. The magnetic response 



of long cylinders ||66|j3^ and narrow rings (the two nontrivial generalizations of one- 
dimensional rings) can be calculated by neglecting the curvature of the circle and solving 
the Schrodinger equation for a rectangle with periodic boundary conditions. Our semiclas- 
sical derivation provides an intuitive and unifying approach to the magnetic response of 
circular billiards and rings of any thickness (for individual systems as well as ensembles) 
and establishes the range of validity of previous studies. Moreover, we present it for com- 
pleteness since it provides a pedagogical introduction to the more complicated ("generic") 
cases of the following sections. 



A. Oscillating density of states for weak field 

By definition, a classical Hamiltonian 7^(p, q) is integrable if there exist as many con- 
stants of motion in involution as degrees of freedom. For bounded systems, this implies (see 
e.g. [37]) that all trajectories are trapped on torus-like manifolds (invariant tori), each of 
which can be labeled by the action integrals 

/^ = ^£pdq (^ = 1,2), (4.1) 

taken along two independent paths Ci and C2 on the torus. (We are dealing with two degrees 
of freedom.) It is moreover possible to perform a canonical transformation from the original 
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(p, q) variables to the action-angle variables (1,0) where I = (/i,/2) and = {ipi^ip2) with 
(^1, in [0, 27r]. Because both, Ji and J2, are constants of motion, the Hamiltonian I2) 
expressed in action-angle variables depends only on the actions. 

For a given torus we note z/j = dH/dIi {i = 1, 2) the angular frequencies, and a = 
the rotation number. A torus is said to be "resonant" when its rotation number is rational 
{a = Ui/u2 where ui and U2 are coprime integers). In that case all the orbits on the torus 
are periodic, and the torus itself constitutes a one-parameter family of periodic orbits, each 
member of the family having the same period and action. The families of periodic orbits can 
be labeled by the two integers (Mi,M2) = {rui,ru2) where (mi,M2) specifies the primitive 
orbits and r is the number of repetitions. Mj {i = 1,2) is thus the winding number of ipi 
before the orbits close themselves. The pair M = (Mi, M2) has been coined the "topology" 
of the orbits by Berry and Tabor. 



For two-dimensional systems, the Berry- Tabor formula can be cast in the form ||5^ , |53 



with 



= ^ dMAE) , (4.2) 

M^{0,0),e 



dmAE) = — — — .. 1/2 ^^^y-J^ ^M- + 7^)- (4.3) 



The sum in Eq. ( [4.2| ) runs over all families of closed orbits at energy E, labeled by their 
topology M (in the first quadrant, that is Mi and M2 are positive integers), and, except for 
self-retracing orbits, by an additional index e specifying tori related to each other through 
time reversal symmetry and therefore having the same topology, gs represents the spin 
degeneracy factor, while Sm.e and tm are, respectively, the action integral and the period 
along the periodic trajectories of the family M. t/m is the Maslov index which counts 
the number of caustics of the invariant torus encountered by the trajectories. For billiard 
systems with Dirichlet boundary conditions, we will also take into account in t/m the phase 
TT acquired at each bounce of the trajectory on the hard walls (and still refer to 57m as the 
Maslov index, although slightly improperly). The energy surface E in action space whose 
imphcit form is ?-^(Ji, J2) = E, is exphcitly defined by the function I2 = gE^h)- We note 
l'^ = (/f^, I^) the action variables of the torus where the periodic orbits of topology M 
live. They are determined by the resonant-torus condition 



a = - 



dli 



^ , (4.4) 



where the first equality arises from the differentiation of T-C{Ii, gE{Ii)) = E with respect to 
Ji. Finally, the last contribution to the phase is given by 7 = sgn{g'^{I^)) . 

The [first] derivation of the Berry- Tabor trace formula follows very similar lines as 
the treatment of the density of states performed in the introduction for the macroscopic 
Landau susceptibility. The EBK (Einstein, Brillouin, Keller) quantization condition is used 
instead of the exact form ( |1.5|) of the Landau levels, followed by the apphcation of the 
Poisson summation rule. While in the latter case this procedure leads to the exact sum of 
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Eq. (|1.6| ), the Berry- Tabor formula is obtained (similar to the treatment of de Haas - van 
Alphen oscillations for a non-spherical Fermi surface) after a stationary-phase approximation 
valid in the semiclassical limit where S ^ h (with a stationary-phase condition according 
to Eq. (0)). 

Given a two-dimensional electron system whose classical Hamiltonian 

^(P,q) = 7^ fp - -A(q)y + l^(q) (4.5) 




2m \ c 

remains integrable for finite values of the transverse field Hz = V x A, the magnetic response 
can be obtained, in principle, from the calculation of the various quantities involved in the 
Berry- Tabor formula at finite fields. However, for weak fields, one can use the fact that the 
field dependence of each contribution cIm to the oscillating part of the density of states is 
essentially due to the modification of the classical action, since this latter is multiplied by 
the large factor 1/h, while the field dependence of the periods and the curvatures of the 
energy manifold can be neglected. Therefore, in this regime we will use for tm and qe the 
values and g% at zero field and consider the first order correction 6S to the unperturbed 
action S^. A general result in classical mechanics |56,^ states that the change [at constant 
energy) in the action integral along a closed orbit under the effect of a parameter A of the 
Hamiltonian is given by 

where the integral is taken along the unperturbed trajectory. Therefore, if the Hamiltonian 
has the form of Eq. ( [4.5|) , classical perturbation theory yields for small magnetic fields H, 

6S = - HA, , (4.7) 

c 

where A, is the directed area enclosed by the unperturbed orbit. This expansion is valid for 
magnetic fields low enough, or energies high enough, such that the cyclotron radius of the 
electrons is much larger than the typical size of the structure (r^ = mcv/eH 3> a, which is, 
e.g., the case for electrons at the Fermi energy in the experiments of Refs. p5| , p6| ). In this 
case we neglect the change in the classical dynamics and consider the effect of the applied 
field only through the change of the action integral. 

For a generic integrable system there is no reason, a priori, that all the orbits of a given 
family M should enclose the same area. However, as pointed out above, a characteristic 
feature of integrable systems is that the action is a constant for all the periodic orbits of a 
given resonant torus. Therefore, the fact that a system remains integrable under the effect of 



a constant magnetic field implies (because of Eq. ([4.7|) ) that all the orbits of a family enclose 
the same absolute area A-M.,f Moreover, since the system is time-reversal invariant at zero 
field, each closed orbit (M, e) enclosing an area ^M,e is associated with a time-reversed 
partner having exactly the same characteristics except for an opposite enclosed area (if the 
orbit is its own time reversal, = 0). Grouping time-reversal trajectories in Eq. ( |4.2| ) at 
if = we have 
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{'^M,e(-^) self-retracing orbits 

£ dli^eiE) = 2 dM,e(^) for non self-retracing orbits • (4-8) 

e=±l 

For weak fields the contribution of self-retracing orbits is unaltered and therefore they do 
not contribute to the magnetic response. For the non self-retracing ones we have 

dm{E,H) = J2 dm,e{E,H) = dl,{E) cos (^^m) , .4m = |^m,.| • (4.9) 
This is the basic relation to be used in the examples that follow. 



B. Circular billiards 



We now apply the preceding considerations to a two-dimensional gas of electrons moving 
in a circular billiard of radius a (where the potential V^(q) is zero in the region |q| < a and 
infinite outside it). Thus we deal with vanishing wavefunctions at the boundary (Dirichlet 
boundary condition). 

In billiards without magnetic field the magnitude p of the momentum is conserved, and 
it is convenient to introduce the wave number, 



since at iJ = the time-of-fiight and the action-integral of a given trajectory can be simply 
expressed in terms of its length L as 

= - L , ^ = kL. (4.11) 
p n 



Following Keller and Rubinow [Q, we calculate the action integrals I = (Ji, I2) by using 
the independent paths Ci and C2 displayed in Fig. |l3(a). The function qe is given by (see 
and Appendix y) 



1 



TT 
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gEih) = - (pay -If - h arccos — , (4.12) 



pa 



where Ji is interpreted as the angular momentum and bounded by < Ji < pa. 

The periodic orbits of the circular billiard are labeled by the topology M = (Mi,M2), 
where Mi is the number of turns around the circle until coming to the initial point after 
M2 bounces. (Obviously M2 > 2Mi.) Elementary geometry yields for the length of the 
topology-M trajectories 

Lm = 2M2a sin (5 , (4.13) 
where 6 = ttMi/M2. The resonant-torus condition, Eq. ( [4.4|) , allows us to obtain as 

/f= pa cos 5, (4.14a) 
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/M = P^^sm5 - 5 cos 5} . (4.14b) 

TT 

The Maslov index of the topology-M trajectories is ffM = 3M2 (M2 bounces, each of them 
giving a dephasing of vr, and M2 encounters with the caustic per period). We therefore have 
all the ingredients necessary to calculate the oscillating part of the density of states at zero 
field: For the non self-retracing trajectories we obtain 

4.(E) ^ /I ^ co= (.Lm + J - |m.) . (4.15) 

The contribution of a self-retracing orbit is just one half of the contribution ( [4.15|) . Its field 



dependent counterpart is obtained from Eq. ( [4.9|) with the area enclosed by the periodic 
orbits given by 

AM = ^^sm26 . (4.16) 

The bouncing-ball trajectories M2 = 2Mi (with zero angular momentum) are self-retracing 
and have no enclosed area; thus they do not contribute to the low field susceptibility. 

Using Eqs. (^]^) and (I^TT^ ), and noting kp = k{ji) = (2/a) {N{fx) /gsY^"^ the Fermi 



wave vector, we obtain the contribution to the magnetic susceptibility associated with AF^^^: 
{kpaf^ X (4.17) 



XL V27r 

,2\2 



Mi,M2>2Mi 



Since we are working with billiards, the temperature factor Rt is given in terms of the 
trajectory length Lm by Eq. (|A5D and the characteristic cut-off length Lc = hvpP/iT. For 
M2 ^ Ml we have Lm — 27iMia and — vrMio^, independent of M2. Performing 
the summation over the index M2 (for fixed value of Mi) by taking the length and area 
dependent terms outside the sum we are left with a rapidly convergent series (whose general 
term is (— 1)^^^/M|). We can therefore truncate the series after the first few terms. In 
Fig. 1^ the sum ( [4.17| ) is evaluated numerically at zero field (solid line) for a cut-off length 
Lc = 6a which selects only the first (Mi = 1) harmonic, and the beating between the first 
few periodic orbits is obtained as a function of wave- vector kp. With only the first two 
primitive orbits (M2 = 3 and 4, dashed line) we give a good account of x^^-* for most of 
the fc-interval. Taking the first four primitive orbits suffices to reproduce the whole sum. 
The short period in kp corresponds approximately to the circle perimeter L = 27ra. Going 
to lower temperatures gives an overall increase of the susceptibility but does not modify 
the structure of the first harmonic contribution since the length of the whispering-gallery 
trajectories is bounded by L. However, for larger values of Lc higher harmonics, namely up 
to Ml of the order of Lc/27ra, will be observed. The predominance of the first few trajectories 
also appears in the beating as a function of magnetic field (not shown) that results from the 
evaluation of (|4.17|) at finite fields. 
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/^From Fig. |^ we see that the susceptibihty of a circular bilhard oscillates as a function 
of the number of electrons (or kp) taking paramagnetic and diamagnetic values. Its overall 
magnitude is much larger than the two-dimensional Landau susceptibility and grows as 
(fcj^a)^/^. We will later show (Sec. |V|) that this finite-size increase with respect to the bulk 
value is distinctive of systems that are integrable at zero field. In order to characterize the 
typical value of the magnetic susceptibility we define 



1/2 



(X«)2 (4.18) 



where, as in section ||, the average is over a kpa interval classically negligible {A{kpa) <^ 
kptt) but quantum mechanically large (A(/c^a) ^ 2n), so that off-diagonal terms 
cos(/c^Lm) cos(A;^Lm') with M 7^ M' vanish under averaging. A remark is in order here 
because at fixed Mi, Lm goes to 2nMia as M2 goes to 00, and {L(^Mi,M2) ~ L(Mi,m^)) can 
be made arbitrarily small by increasing M2 and M2. Therefore, for any interval of kpa over 
which the average is taken, some non-diagonal terms should remain unaffected. Neverthe- 
less, because of the rapid decay of the contribution with M2, these non-diagonal terms can 
be neglected in practice for the experimentally relevant temperatures. The typical zero-field 
susceptibility of the circular billiard is then given by 
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(4.19) 



2 Mi,M2>2A/i -^m/o M| 

Numerical evaluation of the first harmonic (Mi = 1) from (|4.18| ) on the kpa interval of Fig. |^ 
with Lc = 6a gives 2.20(kpay^'^XL (dotted horizontal line), while Eq. ( [4.19| ) restricted to 
M2 < 6 yields 2.16{kpaY^'^XL illustrating the smallness of the off-diagonal and large- M2 
terms. 

For an ensemble made of circular billiards with a dispersion in size or in the number of 
electrons such that A(kpa) > 2ti, the term x^^'^ yields a vanishing contribution to the average 
susceptibility. In such a case it is necessary to go to the next-order free-energy term AF^'^\ 
whose associated contribution x^^^ yields the average susceptibility by means of Eqs. ( p.2|) 
and ( p.20c| ). For the same reason as above one can show that only diagonal terms of [N^^'^Y 
survive the kpa average, in spite of the degeneracy of the length of the closed orbits as M2 
goes to 00. One therefore has 

- kpa 2^ cos \^-^Am j Rt{Lm) ■ (4.20) 



XL VT 



Mi,M2>2Mi 



Again, the terms generally decay rapidly with M2 (as 1/M|), and for a cutoff length Lc 
selecting only the terms with Mi = 1 the total amplitude at zero field (5.2/c^a) can be ob- 
tained from the first few lowest terms. The low field susceptibility of an ensemble of circular 
billiards is paramagnetic and increases linearly with kpa. As for the x^^^ contribution, we 
will show in the sequel that this behavior does not necessitate the integrability at finite 
fields, but rests only upon the integrability at zero field. 

Up to now there have not been measurements of the magnetic response of electrons in 
circular billiards (individual or ensembles). Our typical (Eq. ( |4.19| )) or average (Eq. ( |4.2CI| )) 
susceptibilities exhibit a large enhancement with respect to the bulk values (by powers of 
kpa). Thus it should be possible to detect experimentally these finite-size effects. 
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C. Rings 



The magnetic response of small rings can be calculated along the same lines as in the 
case of the circles. The ring geometry deserves special interest since it is the preferred 
configuration for persistent current measurements. In a ring geometry at = we have two 
types of periodic orbits: those which do not touch the inner disk (type-I), and those which 
do hit it (type-II). (See Fig. |l^ of Appendix 0; we note by a and b respectively the outer 
and inner radius of the ring.) The function QEih) has two branches corresponding to the 
interval to which the angular momentum Ji belongs. For pb < Ii < pa, (type-I trajectories) 
qe has the same form (|4.12| ) as for the circle, while for < Ji < pb, (type-II trajectories) we 
show in Appendix |C| that 
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(4.21) 



The type-I trajectories are labeled in the same way as for the circle by the topology 
M = (Ml, M2) representing the number of turns Mi around the inner circle until returning 
to the initial point after M2 bounces on the outer circle. We therefore obtain the resonant- 
tori condition Eqs. (|4.14| ) and the same contribution ( [4.15| ) to the oscillating part of the 
density of states as in the case of the circle. The only difference is that in the Berry- Tabor 
trace formula (Eq. ( |4.2| )) the sum corresponding to type-I trajectories is now restricted to 
M2 > M2(Mi) = Int[Mi7r/ arccos r]. We note by Int the integer-part function and r = b/a. 
We stress the fact that the minimum value of M2 is itself a function of Mi. The previous 
restriction can also be expressed as cos 5 > r, with 5 = ttMi/M2. Type-II trajectories can 
be labeled by the topology M = (Mi,M2), where Mi is the number of turns around the 
inner circle in coming to the initial point after M2 bounces on the outer circle. We have the 
same restriction M2 > M2(Mi) as for type-I trajectories, and we can use 57m = since there 
are 2M2 bounces with the hard walls and no encounters with the caustic. From ( ^4.21|) we 
obtain the resonant-torus condition 



Jf = pb 
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2r cos 5 



(4.22a) 
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(4.22b) 



The H = contribution to the oscillating part of the density of states from non self- 
retracing type-II trajectories with topology M is given by 



= M 



while its length is 



2 gsO^m [(1 — r cos(5)(r cos5 — r^)] 
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(4.23) 



Lm = 2M2a Vl + r2 - 2r cos 5 . 



(4.24) 
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The small field dependence follows from Eq. ( [4.9|) using the enclosed area 



A 



M 



Moab sin 6 



(4.25) 



In the case of annular geometries it is customary to characterize the magnetic moment 
Ai of the ring by the persistent current 




(4.26) 



T,N 



In order to pass from the applied magnetic field H to the flux $ we use the area A of the outer 
circle ($ = AH, A = t^o?) as defining area. (For thin rings, all periodic orbits with the same 
repetition number Mi enclose approximately the same flux Mi$.) Applying Eqs. ( p^.lScj )- 
( p.20| ), and calling Jq = evpjlT^a the typical current of one-dimensional electrons at the 
Fermi energy, the persistent current of a ring billiard can be expressed as the sum of two 
contributions corresponding to both types of trajectories: 
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(4.27) 
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In Fig. 1^ we present the first harmonic Ii^^ of the persistent current for a thin ring and 
a cut-off length Lc = 6a (solid line). (I.e., we are considering the winding number Mi = 1.) 
The contribution of type-I trajectories (dashed line) is similar as in the case of the circle: 
a rapidly convergent sum showing as a function of kp the beating between the first two 
trajectories (M2 and M2 + 1). On the other hand, Eq. ( [4.28b| ) shows that the trajectories 
with low values of M2 (i.e. M2 ~ M2) contributing to T^jj have negligible weight due to the 
small stability prefactor caused by the defocusing effect exerted by the inner disk (cos 5 ~ r). 
The sum is dominated by trajectories with M2 > M2 and therefore we loose the previous 
beating structure in the total The short period in kp still corresponds to the circle 

perimeter L. 

As in the previous subsection, we characterize the typical value of the magnetic response 
by averaging (J*^^^)^ over a fc^a-interval containing many oscillations, but yet negligible on 
the classical scale. 
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In the same way as for the circular bilhard, one can in practice consider that, despite the 
degeneracy in the length of type-I trajectories for large M2, only diagonal terms (in both 
index M and trajectory-type) survive the averaging for large enough A{kpa). Therefore 
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^ jj) are obtained from Eqs. ( |4.28|) simply by replacing the average 



where (X^,/)^ and (j£^ 
of cos^{kpL^ + 7r/4 - 3M27r/2) and sin^(fc^LM + 7r/4) by 1/2. 

In Fig. ^ we present the typical persistent current and its two contributions for various 
ratios r = h/a and cut-off lengths for the first harmonic (Mi = 1). The contribution 
^M7 of type-I trajectories dominates for small r (where the inner circle is not important 
and we recover the magnetic response of the circular billiard) while type-II trajectories take 
over for narrow rings. The crossover r depends on temperature through due to the 
different dependence of the trajectory length on M (Eqs. ( |4.13 ) and ( 4.24| )) for both types 
of trajectories. 

As in the case of x^^^ for the circular billiard, J^^-' gives a vanishing contribution to the 
persistent current of an ensemble of rings with different sizes or electron fillings as soon as 
the dispersion in k^a is of the order of 27r. We therefore need to go to the term AF*^^^ in 
the free-energy expansion, which is obtained (see Eq. ( 2.20c| )) from 
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where N^jifi) and NMjiifi) are given in terms of the respective contributions to the field 
dependent density of states through Eq. ( p.l5b|) . For an ensemble with a large dispersion of 
sizes only diagonal terms survive the average and we have (with D = gs'mA(l — r^)/(27r/i )) 
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The kp dependence of the average persistent current is linear (through Jq), similarly to 
the case of the average susceptibility of an ensemble of circular billiards. 
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Thin rings 



In the case of thin rings (a ^ 6, r ~ 1) further approximations can be performed 
on Eqs. ( |4.28|) and ( |4.33| ) using (1 — r) as a small parameter, giving more compact and 
meaningful expressions for the typical and average persistent currents. Since in addition 



this is the configuration used in the experiment of Ref. p6[, we shall consider more closely 
this limiting case. First, we note that 5 



arccos r 



2 1 - r < 1. Thus 
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and for M2 > M2, the area and length of contributing orbits can be approximated by 

^ ^ MxA = Mina^ ; Lm ^ MiL = Mi2na . (4.35) 

For the length of type-II trajectories we have Lm — MiL for M2 ~ M2, and Lm ~ 2M2(a— 6) 
when M2 ^ M2. AH trajectories with winding number Mi enclose approximately the same 
flux Mi$, and the field dependent terms in Eq. (|4.27| ) may be replaced by sin (27rMi$/$o)- 
There is therefore no difference between the case that we study (where a uniform magnetic 
field H is applied) and the ideal case of a flux line $ through the inner circle of the ring. 
The length dependent factors can also be taken outside the sum over M2 since the main 
contribution of type-II trajectories comes from M2 ~ 7rMi/[5^/^(l — r)^/^]. Even if these 
M2's are much larger than M2, their associated Lm are still of the order of MiL to leading 
order in 1 — r. 

Turning now to the typical and ensemble average currents, it should be stressed that 
for narrow rings it is necessary to go to fairly large energies before an average on a scale 
being quantum mechanically large but classically small is possible. Indeed, one has for both 



types of trajectories kp{Lj^^^^ — L 
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{AV2/3)7rAfy/T^, where 

A/" = kp{a — b)/TT is the number of transverse occupied channels. Therefore, Af should be 
much larger than (1 — r)^^/"^ if one wants to assume A{kpa) sufficiently large to average 
out all non-diagonal terms without violating the condition A{kpa) <^ kpa. Supposing the 
previous condition is met, and introducing the typical amplitudes jT^^j j and jT^^j jj of each 
harmonic, we write 
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Since M2 ^ 1 we can convert the previous sums into integrals and obtain 
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In leading order in 1 — r the persistent current is dominated by type-II trajectories 
(independent of the temperature) and given by 



^ = ^ gs VAT 
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(4.39) 



consistent with the result of Ref . . For the next order term the contribution from type-I 
trajectories is cancelled by that of type-II resulting in the relatively flat character of the 
curves for J*^*) in Fig. |. 

For the current of an ensemble of thin rings, the calculations are similar to those of 
Eqs. (|4.38| ), and in leading order in 1 — r we obtain: 



(4.40) 
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Type-II trajectories once again dominate the average magnetic response of thin rings and 

the amplitude for the first harmonic is if^ / Iq — (2gs/7r^) sin {Att^ /^q)E^{L), independently 
of the number of transverse channels M . The average persistent current shows the halfing 
of the flux period with respect to J*-^^ characteristic for ensemble results (as found in the 
disordered case and consistently with the results for averages in the following sections). 



Comparison with Experiment 

Persistent currents have been measured by Mailly, Chapelier and Benoit in a thin 
semiconductor ring (with effective outer and inner radii a = 1.43/im and b = 1.27/im) in 
the ballistic and phase-coherent regime (/ = 11/im and L$ = 25/im). The Fermi velocity is 
Vp = 2.6 X lO^cm/s and therefore the number of occupied channels is A/" — 4. The quoted 
temperature of T = IbmK makes the temperature factor irrelevant for the first harmonic 
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{Lc — 30a, Rt{L) ~ 1). The magnetic response exhibits an hc/e flux periodicity and changes 
from diamagnetic to paramagnetic by changing the microscopic conflguration, consistently 
with Eqs. ( [4.27| )-( P:.28| ). Unfortunately, the sensitivity is not high enough in order to test the 
signal averaging with these microscopic changes. The typical persistent current was found 



to be An A, while Eq. ( ^.391 ) and Ref. [^] would yield In A. The difference between the 
theoretical and measured values is not signiflcant given the experimental uncertainties as 
discussed in Refs. and P^. Moreover, as we stressed above, a very large kpa interval is 
needed for the average of (/W)2 order to recover J*^*); otherwise we expect large statistical 
fluctuations. As in the case of the susceptibility of squares that we analyze in the next 
section, residual disorder (reducing the magnetic response without altering the physical 
picture) and interactions may be necessary in order to attempt a detailed comparison with 
the experiment. Clearly, new experiments on individual rings of various thickness and on 
ensembles of ballistic rings would be helpful in order to test the ideas of the present section. 
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V. SIMPLE REGULAR GEOMETRIES: THE SQUARE 



The circular and annular billiards studied in section ^ have the remarkable property 
that, due to their rotational symmetry, they remain integrable under the application of a 
magnetic field. However, for a generic integrable system (a regular geometry) any pertur- 
bation breaks the integrability of the dynamics. Moreover, the periodic orbits which are 
playing the central role in the semiclassical trace formulas are most strongly affected by the 
perturbation. Indeed, the Poincare-Birkhoff theorem |^ states that as soon as the magnetic 
field is turned on, all resonant tori (i.e. all families of periodic orbits) are instantaneously 
broken, leaving only two isolated periodic orbits (one stable and one unstable). It is there- 
fore no longer possible to use the Berry- Tabor semiclassical trace formula to calculate the 
oscillating part of the density of states for finite field since it is based on a sum over resonant 
tori, which do not exist any further. One has therefore to devise a semiclassical technique 
allowing to calculate d°^'^{E) for nearly, but not completely, integrable systems. 

To achieve this, it is necessary to go back to the basic equations from which the standard 
semiclassical trace formulae (Gutzwiller Balian-Bloch ||5^, Berry- Tabor [^) are de- 
rived. The density of states d{E), Eq. (pT^), is related to the trace of the energy dependent 
Green function G{q, q'; E) by 
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(5.1) 



where again the factor gs = 2 comes from the spin degeneracy. G(q, q'; E) has a singularity 
(logarithmic in two dimensions) when r ^ r' which just gives the smooth [Weyl] part d{E) 
of the density of states in a leading order semiclassical expansion. However, in order to 
consider only the oscillating part of d{E) one can use the semiclassical approximation of the 
Green function 
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where the sum runs over all classical trajectories t joining q and q' at energy E. St is the 
action along the trajectory t, Dt a determinant involving second derivatives of the action 
(the general expression of which is given in appendix |^) and rjt is the Maslov index of the 
trajectory, i.e. the number of focal points encountered when traveling from q to q'. As in 
section we shall also take into account in rjt the phase vr acquired at each refiection at 
the wall of a billiard with Dirichlet boundary conditions. 

By taking the trace ( |5.1|) the sum in Eq. ( p.2|) becomes a sum over all orbits closed in 
configuration (i.e. q) space, to which we will refer in the following as recurrent orbits. The 
standard route to obtain d°^'^ is to evaluate this integral by stationary-phase approximation. 
This selects the trajectories which are not only closed in configuration space (r' = r), but 
also closed in phase space (p' = p), i.e. periodic orbits. When these latter are [well] isolated 
the Gutzwiller Trace Formula is obtained. For integrable systems, all recurrent orbits 
are in fact periodic since the action variables are constants of motion. 

Periodic orbits appear in continuous families associated with resonant tori. All orbits 
of a family have the same action and period, and one can calculate the density of states 
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using the Berry- Tabor Formula as described in the previous section. For systems such as the 
square bilhard, the physical effect which generates the susceptibility comes along with the 
breaking of the rational tori, so that just ignoring this, i.e. using the Berry- Tabor Formula, is 
certainly inadequate. On the other hand, for — the remaining orbits are not sufficiently 
well isolated to apply the Gutzwiller Trace Formula. Therefore, as stated before, we need a 
uniform treatment of the perturbing field, in which not only the orbits being closed in phase 
space are taken in account, but also the orbits closed in configuration space which can be 
traced back to a periodic orbit when H 0. 

In this section we show how this can be performed in the particular case of a square bil- 
liard. Because of the simplicity of its geometry, the integrals involved in the trace Eq. ( p.l|) 
can be performed exactly for weak magnetic fields. Moreover, the square geometry deserves 
special interest since it was the first microstructure experimentally realized to measure the 
magnetic response in the ballistic regime. We present here a semiclassical approach address- 
ing the physical explanation of the experimental findings of Ref. which have pointed 
the way for the ongoing research. In order to obtain semiclassical expressions for the suscep- 
tibility of individual and ensembles of squares we will proceed as outlined in section ^ We 
will calculate the density of states and use the decomposition of the susceptibility according 
to Eq. ( ^.191 ) into contributions corresponding to AF^^^ and AF^'^\ In section |V| we present 
the theory for a generic integrable system perturbed by a magnetic field, generalizing the 
results of this section. 



A. Oscillating density of states for weak field 

To start with, we consider a square billiard (of side a) in the absence of a field. Each 
family of periodic orbits can be labeled by the topology M = {M^, My) where and My 
are the number of bounces occurring on the bottom and left side of the billiard (see Fig. 
The length of the periodic orbits for all members of a family is 



2a JM2 + M2 . (5.3) 



The unperturbed action along the trajectory is, as for any billiard system, S^/h = kLf^ 
where k is the wavenumber. The Maslov indices are ?7m = ^{M^+My), and we will omit them 
from now on since they only yield a dephasing of a multiple of 27r. Finally the unperturbed 
determinant reduces to 

Tfl 

One way to obtain this result is to use the method of images (see Fig. ^) and express the 
exact Green function (^(q, q'; E) in terms of the free Green function (^^(q, q'; E) as 



G(q, q'; E) = G°(q, q'; + ^ e,G°(q„ q') , (5.5) 



where the q^ represent all the mirror images of q by any combination of symmetry across 
a side of the square, and = +1 or — 1 depending on whether one needs an even or odd 
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number of symmetries to map q on q^. G"^(q, q'; E) gives the above mentioned logarithmic 
singularity of G when q' — > q, but the long range asymptotic behavior of the two-dimensional 
free Green function 

G^(a. a') 1 ^ exp(zfc|q--q,|) 



can be used for all other terms (images). 

For sufficiently weak magnetic fields, one may follow the same approach as in the previous 
section, keeping in Eq. ( |5.2| ) the zero'th order approximation for the prefactor -Dm, and using 
the first-order correction 6S to the action which, as expressed by Eq. ([4.7|) , is proportional 
to the area enclosed by the unperturbed trajectory. Here however, as is the generic case 
(and contrary to circular or annular geometries) the area enclosed by an orbit varies within 
a family. 

Let us consider the contribution to the density of states of the family of recurrent trajec- 
tories which for H^O tends to the family of shortest periodic orbits with non-zero enclosed 
area, that plays a crucial role in determining the magnetic response, as already recognized 
in Ref. p5[. For H = 0, this family consists in the set of orbits which, say, start with an 



angle of 45 degrees with respect to the boundary on the bottom side of the billiard at a 
distance xq (0 < < a) from its left corner, bounce once on each side of the square before 
returning to their initial position (family M = (1, 1), see Fig. §(a)). It is convenient to use 
as configuration space coordinates xq which labels the trajectory, the distance s along the 
trajectory, and the index e = ±1 which specifies the direction in which the trajectory is 
traversed. In this way, each point q is counted four times corresponding to the four sheets 
of the invariant torus. The enclosed area Ae{xQ,s) obviously does not depend on s and is 
given by 

•Aei^o) = e 2 Xq {a — Xq) . (5.7) 

Periodic orbits are those paths for which the action is extremal (VS = p' — p = 0). Therefore 
Eqs. ([4.7|) and (|5.7| ) illustrate the contents of the Poincare-Birkhoff theorem, that for any 
non-zero field only the two trajectories corresponding to Xq = a/2 remain periodic (one 
stable, one unstable according to the two possible directions of traversal). The contribution 
of the family (1,1) to (i°'^^(E) is du{E) = -(gs/vr) Im GuiE). Inserting Eqs. (^) and 
into the integral of Eq. (^.1|) we have 




e=±l 



., .2ee , . 

ikLii + t-—Hxo[a — Xq) 
he 



(5.8) 



The contribution to the density of states of the family (1,1) factorizes into an unperturbed 
(Berry- Tabor-like) term and a field dependent factor 

dn{E,H) = dUE)C{H) (5.9) 

with 

* - <in(i/ = 0) = f Jf^^^^. - («n + j) , (5.10) 
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and 

'2e„ , A 1 



\ ra. /2e \ 1 

CiH^) = - dxQ cos I —Hxo{a — xq) ) = [cos{7!'(p)C{y/7np) + sm{TT(p)S{^J7i(p)] . 
a Jo \Tic J v2o9 



2^ 

(5.11) 



C and S respectively denote the cosine and sine Fresnel integrals |70], and 



»> = ^ (5.12) 

is the total flux through the square measured in units of the flux quantum ($o = hc/e). 
For the circular and annular geometries, the fleld dependence of the density of states, and 
therefore the susceptibility, was related to the dephasing between time reversal families of 
orbits. Here, Eq. ( p.ll| ) expresses that the dependence of d°^^ on the fleld is also determined 
by the fleld induced decoherence of different orbits within a given family. 

As soon as (p reaches a value close to one, the Fresnel integrals can be replaced by their 
asymptotic value 1/2, which amounts to evaluate C{(p) by stationary phase, i.e. 

^c, , cos(7ro9 — 7r/4) 

C\v) = r^' ■ (5.13) 

This means that for ip > 1 the dominant contribution to C{ip) comes from the neighborhood 
of the two surviving periodic orbits {xq = a/2,e = ±1), and the oscillations of C{ip) are 
related to the successive dephasing and rephasing of these orbits. In fact, one would have 
obtained just df^ = dl-^^C^^ip) by evaluating the contribution to the density of states of the 
two surviving periodic orbits using the Gutzwiller trace formula with a flrst-order classical 
perturbative evaluation of the actions and stability matrices. C^{f) however diverges when 
while the full expression Eq. ( p.llD simply gives C(0) = 1. 
To compute the contribution dm of longer trajectories, it is worthwhile to write {M^, My) 
as {rux, ruy), where Ux and Uy are coprime integers labeling the primitive orbits and r is the 
number of repetitions. As illustrated in Fig. ^(b), for any orbit of the family the square can 
be decomposed into Ux x Uy cells, such that the algebraic area enclosed by the trajectory 
inside two adjacent cells exactly compensate. Therefore, keeping xq as a label of the orbit 
(with Xq G [0,a/ux] to avoid double counting), the total area enclosed by the trajectory 
{rux,ruy) is 

{0 Ux or Uy even 

r^^^ ..and ..odd ' ^'''^^^ 

where Ae{xo) is given by Eq. (|5.7| ). From the above equation, and proceeding in the same 
way as for the orbit (1,1) Eq. ( p.9| ) can be generalized to 

d-mi^) '^x or Uy even 

d^{E^ H) = { ^0^^^^ ^(np\ , (5.15) 

yUxUy J 
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where C{(p) is given by Eq. (|5TT|) and c?^ = dMiH = 0) is the zero-field contribution of the 
family M 



B. The susceptibility: individual samples vs. ensemble averages 

For clarity of the presentation we will calculate in a first stage the susceptibility contri- 
bution of the family (1,1) of the shortest fiux enclosing orbits only. This corresponds to the 



temperature regime of the experiment Ref. where the characteristic length Lc given by 



Eq. ( |X5| ) is of the order of Ln, the length of the shortest orbits, and contributions of all 
longer orbits are eliminated due to temperature damping. In the next subsection we will 
state the results valid at arbitrary temperature by taking into account the contribution of 
longer orbits. 

iFiom. the expressions (^.9|) and ( ^.10|) of the contributions of the family (1, 1) to 



d"^'^{E,H) one obtains the corresponding contribution to AF*^^^ (Eqs. ( |2.15c| ) and ( |2.20b| )) 



as 



AFfl\H) = ^ ( ) {k,af^'' sin ( k,Ln + ^ ) C{H)Rt{U^) . (5.17) 




TT 



Rt{Lii) is the temperature dependent reduction factor Eq. (|A5|) , valid for billiard systems. 
The field-dependent factor C{(p) is given by Eq. ( |5.11| ). Taking the derivatives with respect 
to the magnetic field, we have [for Lc — Lu] 



The susceptibility of a given square oscillates as a function of the Fermi energy and can 
be paramagnetic or diamagnetic (see Fig. §|(a)). Since we are considering only one kind of 
trajectory the typical susceptibility x*^*"* (with the definition ([4.18|) ) is simply proportional 
to the prefactor of x^^^- Therefore, it is of the order of (kpa)^^'^, which is much larger than 
the Landau susceptibility Xl- As shown in Fig. ||b (solid line) x^^^ exhibits also (by means of 
d'^C/d(f'^) oscillations as a function of the fiux at a given number of electrons in the square. 
The divergent susceptibility obtained from (dashed line) provides a good description of 
X^-*^) for (f ^ 1. 

For a measurement made on an ensemble of squares of different sizes a, x^^^ vanishes 
under averaging if the dispersion of kpLn across the ensemble is larger than 27r. In that 
case the average susceptibility is given by the contribution to AF^^^ arising from the (1, 1) 
family (Eq. ( |2.20c| )). Proceeding in a similar way as for the first-order term, the contribution 
of the family (1, 1) to the integrated density A^"^'^ is given by Eq. ( |2.15b| ) as 




Num,H) = -gs ( ) (k^aY^-' cos ( kpLu + ^ ) C{H)Rt{L^^) . (5.19) 
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To calculate x^^^ we have to consider AF(2) = {N^^^f /2D (with D = {gsma'^) / {27ih'^)) , and 
in particular the term 

{Nuifi, H)f g^h' ^^^2 f ^^^^^ ^ ^\ ^2(^) ^2 _ (5 20) 



2D (V2)37r2ma2 ^ V 4 

This contribution is of lower order in kpa than that of AFil\ but its sign does not change as 
a function of the phase kpLn. Therefore the squared cosine survives the ensemble average]^ 
and we obtain, performing the derivatives with respect to ip (still in the regime Lc — Lu), 



^ kpa-—RULn). (5.21) 



Xl (v^7r)3 
The total averaged susceptibility is therefore 



X = -Xl + X^^^ 



since, as seen in section III, one has also to include the diamagnetic (bulk) "Landau term" 



—Xl arising from h corrections to In the regime L c:^ we are considering here, Xl 
is negligible with respect to x^^) _^ ^^d one can use x — X^^"*- Note however that 
when Lc -C L, Eqs. ( |5.18| ) and ( [5.21D remain valid but x^^^ ^is well as x*-^^ is exponentially 



suppressed. In this "trivial" regime x (and thus x) reduces to the Landau susceptibility, 
and becomes independent of the underlying classical dynamics. The linear dependence of 
the average susceptibility on hp is shown in Fig. |^. 

Since C has its absolute maximum at if = 0, the average zero-field susceptibility is para- 
magnetic and attains a maximum value of ||57| , |59 



4^2 



X(^KH = 0) = ^ kpa Xl RULu) . (5.22) 

For small fields the average susceptibility (thin solid line. Fig. ^b) has an overall decay as 
l/(p and oscillates in sign on the scale of one flux quantum through the sample. As in 
the disordered case the period of the field oscillations of the average is half of that of 
the individual systems (see Fig. ^(b)). In our case the difference can be traced to the 
dependence that appears in Eq. ( |5.21| ) in contrast to the simple C dependence of Eq. ( p.lSp 



For an ensemble with a wide distribution of lengths (as in Ref. ||2^) an average (■ ■ ■) on 
a classical scale (i.e. Aa/a <^ 1) rather than on a quantum scale (A(fc^a) ~ 2tt) needs to be 
performed, and the dependence of C on a (through ip) has to be considered. Since the scale 



^ Beside the orbits (1,1) the orbits (1,0) and (0,1) which are even shorter contribute to AF^^^ in 
the limit ~ Lu. Since they do not enclose any flux the second derivative of ^F^q^ with respect 
to H, i.e. xio^ can be neglected for small fields. However, they enter into x^^^ by means of the cross 
products (A''io + A''oi)A^ii in {N°^^)'^. Nevertheless, they play no role for the averaged x^^^ because 
A'^io and A^n do not oscillate with the same frequency and therefore their product averages out. 
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of variation of C with a is much slower than that of sin^ {kpLu) we can effectively separate 
the two averages and obtain the total mean by averaging the local mean: 



(X) = dax P{a) , (5.23) 



where the quantum average x is given by Eq. (^.21|) and P{a) is the probability distribution 



of sizes a. Taking for P{a) a Gaussian distribution with a 30% dispersion we obtain the 
thick solid line of Fig. ^b. The low-field oscillations with respect to (f are suppressed under 
the second average, while the zero-field behavior remains unchanged. 

The expected value for the susceptibility measured in an ensemble of n squares is n{x) oc 
nkpa, with a large statistical dispersion of \/nx^^^ oc y/n^kpa)^^"^ . However, for experiments 
like the one of Ref. where n ~ 10^ ^ kpa ~ 10^, it is not possible to obtain a diamagnetic 



response by a statistical fluctuation. 



C. Contribution of longer orbits 

In the zero temperature limitQ or more generally if one is interested in results valid 
at any temperature, it is necessary to take also into account the contribution of longer 
trajectories. This can be done following exactly the same lines as for the contribution of 
the family (1,1). From Eqs. ( ^.15|) and (|5.16| ) one obtains the contribution of the family 

M = {Mx,My) = {rux,ruy), (where Ux and Uy are coprime) to AF^^^ 

.2 / o3„ \V2 



^(i)^m _ I ^ .^3/2 fu r.. , ^ 



AF^\H) = ^ — ^ (M)'/' sin A;^Lm + - ) CM(<^)/?r(^M) , (5.24) 



where 




Ux or Uy even 

Cm(^) „^ ■ (5.25) 



Lm and the function C{ip) are given respectively by Eqs. (|5.3|) and ( ^.11|) . In order to get 
X^^) we have to take the second derivative of Cm with respect to the magnetic field. This 
yields zero if either Ux or Uy is even and a factor r"^ / {uxUyY , if both are odd. We therefore 
obtain 



— = %^ {kpafl^Y. 1/2/ 2 I — IwIT ^^^^^ f^F^M + 7) Rt{L 

Xl TT^/^ r r^l\ui, + Uiifl^{UxUyY \ A J \UxUyJ 



Mi 



odd 

(5.26) 



It should be kept in mind however that the expansion in Eq. (2.19) is a priori not valid when 
T ^ 0. 
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valid at any temperature. 

The low temperature result for x^"^^ follows in essentially the same way, but taking the 
average is made rather intricate in the case of a square (as compared for instance to a 
rectangle) because of the degeneracies in the lengths of the particular orbits of this system. 
Indeed, there are infinitely many integers which can be decomposed in at least two different 
ways into sums of two squares. For instance, 11^ + 7^ = 13^ + 1^ = 170. As a consequence, 
Lii 7 = Li3^i, and NnjNi^^i ^ 0. An explicit formula for x(^) therefore requires to handle 
correctly all the non-diagonal terms containing orbits of degenerated lengths which do not 
average to zero. This leads to a number theoretical problem (i.e. characterizing all numbers 
which decomposition as the sum of two squares is not unique), with which we do not deal and 
which moreover will be seen to be of no practical relevance. Therefore, instead of considering 
a square, we will give the expression for x^^) for a rectangle of area a} and of horizontal and 
vertical lengths a ■ e and a ■ e~^. In that case, all the formulae given in section |V A| remain 
valid. As the only difference one has now 

Lm = 2a^(M,/e)2 + (M,e)2 

instead of Eq. ( ^.3|) , which does not give rise to length degeneracies if, as we will suppose, 
e"^ is irrational. Noting that the prefactor of depends as on the length of the orbit 
(instead of L^^"^ for ^F^), one obtains for the canonical correction to the susceptibility 



^ = -lfc,aEE ^ 372 (CT(^) Rl{L^). (5.27) 

odd 



The equations ( p.26[ ) and ( |5.27| ) show that even at zero temperature the strong flux 



cancelation typical for the square (or rectangular) geometry generates a very small prefactor 
I / {r^/'^{u1 + u^y)^/'^{uxUyY) for x^^^ (square geometry) and l/(r {{ux/e)"^ + {uyc)'^)^^'^ {uxUyY) 
for x*-^^ (rectangular geometry). For the second shortest contributing primitive orbit, M = 
(1,3), this yields for instance for X(^) a damping of 1/(9 x 10^/^) ~ 0.0062. For x^^^ the 
multiplicative factor is even smaller. In practice only the repetitions (r, r) of the family 
(1,1) will contribute significantly to the susceptibility, and one can use Eqs. (|5.26| ) and 



( p.27| ) keeping only the term Ux = Uy = 1 oi the second summation. As a consequence, all 
the complications due to the degeneracies in the length of the orbits for the square are of no 
practical importance (Eq. ( p.27| ) restricted to m^, = = 1 can be used for the square with 



e = 1), showing why their detailed treatment was not necessary. As illustrated in Fig. 10 
for x^'^\ the repetitions of the orbit (1,1) are yielding a diverging susceptibility at zero field 
when the temperature goes to zero, but barely affect the result even as T ^ for finite if, 
where the contributions of the repetitions do no longer add coherently. 



D. Numerical calculations 

As a check of our semiclassical results we calculated quantum mechanically the orbital 
susceptibility of spinless particles in a square potential well [—a/2, a/2] in an homogeneous 
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magnetic field. Within the symmetric gauge A = H{—y/2, x/2, 0) the corresponding Hamil- 
tonian in scaled units x = x/a and E = [ma?' /fi^)E reads 

« = - 2 i» ^ j " " " j ^ ' ^ * ' ■ 

with the normalized flux if defined as in Eq. (|5.12| ). Taking into account the invariance of the 
Hamiltonian ( ^.281) with respect to rotations by vr, 7r/2 we use linear combinations of plane- 



waves which are eigenf unctions of the parity operators P,r, P7r/2, respectively. Omitting the 
tilde in order to simplify the notation, they read 

V2[S^{x)C^{y)±iC^{x)Sn{y)] ; (P. = -1) , (5.29) 
V2[Cn{x)C^{y) ± C^{x)Cn{y)] 



\f2i \Sn[x)Sra{y) ± S.m{x)Sn{y)\ 



(P. = +1) (5.30) 



with Sn{u) = sin(n7ra), n even, and Cm{u) = cos(m7rn), m odd, obeying Dirichlet boundary 
conditions. In this representation the resulting matrix equation is real symmetric and de- 
composes into four blocks representing the different symmetry classes. By diagonalization we 
calculated the first 3000 eigenenergies taking into account up to 2500 basis functions for each 
symmetry class. A typical energy level diagram of the symmetry class {P-n,PTT/2) = (1; 1) 
as a function of the magnetic field is shown in Fig. |l|(a). In between the two separable 
limiting cases = and (p — >■ 00 the spectrum exhibits a complex structure typical for a 
non-integrable system which classical dynamics is at least partly chaotic. 

We calculate numerically the (/rand- car^on^ca/ susceptibility (Eq. ( |1.2|) , Fig. |l|b) from 



^"^'") = --2^2 E T- ur. ^ (5.31) 

a'^oH^ ~^ 1 + exp[p(ei — /i)J 



where gs accounts for the spin degeneracy and ej denotes the single particle energies. 



However, in order to address the experiment of Ref. |^ and to compare with the semi- 
classical approach of the preceeding subsection we have to work in the canonical ensemble. 
At T = the free energy F reduces to the total energy and the canonical susceptibility 
(Eq. ( |2.2| )) is given as the sum 

X^-0)--|||^ (5.32) 

over the curvatures of the N single particle energies e^. The susceptibility is therefore domi- 
nated by large paramagnetic singularities each time the highest occupied state undergoes a 
level crossing with a state of a different symmetry class or a narrow avoided crossing with 
a state of the same symmetry. This makes T = susceptibility spectra of quasi-integrable 
billiards (with nearly exact level crossings) or systems with spectra composed of energy 
levels from different symmetry classes (as it is the case for the square) looking much more 
erratic than those of chaotic systems with stronger level repulsion p7| . 

The T = peaks are compensated once the next higher state at a (quasi) crossing is 
considered, and therefore disappear at finite temperature when the occupation of nearly 
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degenerated states becomes almost the same. Thus finite temperature regularizes the sin- 
gular behavior of x at T = and of course describes the physical situation. We obtain the 
canonical susceptibility at finite T from 

The canonical partition function Zn{P) is given by 

ZN{f3) = E exp[-/3K(iV)] (5.34) 

{a} 

with 

oo oo 

E^iN) = J2eiK , N = J2nt. (5.35) 
i=i i=i 

The n° G {0, 1} describe the occupation of the single particle energy-levels. A direct nu- 
merical computation of the canonical partition function becomes extremely time consuming 
at finite temperature. We approximate the sum in Eq. ( |5.34 ) which runs over all (infinitely 



many) occupation distributions {a} for N electrons by a finite sum Zjq{M] (3) over all pos- 
sibilities to distribute N particles over the first M levels with M > N sufficiently large. 
Following Brack et al. |7TI| we calculate Zn{M; (3) recursively using 



Zn{M- P) = Zn{M -1;(3) + Zn-i{M - 1; (3) exp(-/5eAf) (5.36) 
with initial conditions 

Zo{M;(3) = l , Z^{N-1;(3) = (5.37) 

and increase M until convergence of Zn{M,P), i.e. the difference between Zp^{M;(3) and 
Zn{M — 1] (3) is negligible. This recursive algorithm reduces the number of algebraic oper- 
ations to calculate Z^ drastically and is fast and accurate even if fc^T is of the order of 10 
or 20 times the mean level spacing, i.e., in a regime where a direct calculation of Z^- is not 
feasible. 



E. Comparison between numerical and semiclassical results 

Our numerical results for the susceptibility of individual and ensembles of squares are 
displayed as the dashed lines in Figs. |^ and |^ and are in excellent agreement with the 



semiclassical predictions of Sec. VB . Fig. m shows the numerical result for the canonical 

— (1) 

susceptibility and the semiclassical leading order contribution Xu at zero field as a function 



of kpa {^AirN/gs in terms of the number of electrons). The temperature kgT is equal to five 
times the mean level spacing A of the single particle spectrum. The quantum result oscillates 
with a period tt/a/^ as semiclassically expected (Eq. ( p.l8| )) indicating the dominant effect 



of the fundamental orbits of length Ln = 2\/2a. The semiclassical amplitudes (solid line) 
are slightly smaller than the numerics because only the shortest orbits are included. 
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Fig. shows the flux dependence of x fo^' ^ fixed number of electrons ~ llOOgs. The 
semiclassical prediction (Eq. (|5.18 ), sohd curve) is again in considerable agreement with the 



quantum result while the analytical result (Eq. ( |5.13| ), dashed line) from stationary phase 
integration yields an (unphysical) divergence for as discussed in Sec. [VB . 



For the numerical calculations we can perform the ensemble average directly and we 
obtain the averages on the quantum scale (thin dashed line, Fig. ^o) or classical scale (thick 
dashed line) by taking a Gaussian distribution of sizes with respectively a small or large Aa/a 
dispersion. Fig. ||(a) depicts the kpa dependence of x assuming a Gaussian distribution of 
lengths a with a standard deviation Aa/a ~ 0.1 for each of the three temperatures kgT/A = 
2, 3, 5. The dashed curves are the ensemble averages of the quantum mechanically calculated 
entire canonical susceptibility x- The dotted lines are the exact (numerical) results for the 

averaged term Xqm = (A^°'^'^qm)^/2A. They are nearly indistinguishable (on the scale of the 
figure) from the semic/assica/ approximation of Eq. (|5.21|) (solid line). Although a small fiux 
if ^ 0.15 has been chosen (here the contribution from the next longer orbits (2,2) nearly 
vanishes) the precision of the semiclassical approximation based on the fundamental orbits 
(1,1) is striking. The difference between the results for x and x^^) gives an estimate for the 
precision of the thermodynamic expansion Eq. ( p.l9| ). The convoluted semiclassical result 
has been shifted additionally by —Xl to account for the diamagnetic Landau contribution 
and is again in close agreement with the numerical result of the averaged susceptibility x- 



F. Comparison with the experiment 

In a recent experiment. Levy et al. measured the magnetic response of an ensemble 



of 10^ microscopic billiards of square geometry lithographically defined on a high mobility 
GaAs heteroj unction. The size of the squares is on average a = 4.5/im, but has a large 
variation (estimated between 10 and 30%) between the center and the border of the array. 
Each square can be considered as phase-coherent and ballistic since the phase-coherence 
length and elastic mean free path are estimated, respectively, to be between 15 and 40 jim 
and between 5 and 10 fim. 

Therefore, it is worthwhile to compare the observed magnetic response with the pre- 
diction of our clean model of non-interacting electrons, to see whether this simple picture 
contains the main physical input to understand the experimental observations, although 
one should control in addition that the residual impurities do not alter fundamentally the 



magnetic response of the system. This is the subject of a forthcoming article Ongoing 
calculations including (weak) disorder indeed indicate that the underlying physical picture 
remains correct. 

At a qualitative level, a large paramagnetic peak at zero field has been observed in 
Ref. [^, two orders of magnitude larger than the Landau susceptibility, decreasing on a 



scale of approximately one fiux quantum through each square. Since there is a large disper- 
sion of sizes we do not observe the field oscillations of the quantum average ( [5.21| ), but the 



comparison has to be established with the classical average resuhs Eq. ( ^ ). The corre- 
sponding results from our semiclassical calculations (Eq. ( p.21| , p.23D ) and the full quantum 
calculations are shown in Fig. ||b) as the thick full, respectively dashed, lines (denoted by (x) 
in the figure). The offset in the semiclassical curve with respect to the quantum mechanical 
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curve is due to the Landau susceptibility Xl and additional effects from bouncing-ball orbits 
(see section |VI1] A) not included in the semiclassical trace. Our theoretical results for the 
flux dependence of the average (x) with respect to a wide distribution in the size of the 
squares agree on the whole with the experiment. However, the diamagnetic response for 
(x) that we obtain for ip ^ 0.5 is not observed experimentally, indicating that there may 
be a more important size-dispersion than estimated. As will be discussed in more detail 
in section |V|, a very large distribution of lengths enhances the effect of the breaking of 
time reversal invariance due to the magnetic fleld, yielding a vanishing average response at 
finite field and a paramagnetic susceptibility at zero field decaying on a fleld scale $0 by the 
dephasing of the contribution of time reversal symmetric orbits to the density of states. 

More quantitatively, the experiment of Ref. yielded a paramagnetic susceptibihty 
at i7 = with a value of approximately 100 (with an uncertainty of a factor of 4) in units 
of Xl- The two electron densities considered in the experiment are 10^^ and 3 x lO^^cm"^ 
corresponding to approximately 10^ occupied levels per square. Therefore our semiclassical 
approximation is well justifled. For a temperature of AOmK the factor 4a/2/ {57i)kpaR'^{Lii) 
from Eq. ( |5.22| ) gives zero fleld susceptibility values of 60 and 170, respectively, in reason- 
able agreement with the measurements. In order to attempt a more detailed comparison 
with the measurements we need to incorporate the suppression of the clean susceptibility 
by the residual disorder, which depends on the strength and correlation length of the im- 
purity potential The fleld scale for the decrease of {xiv)) is of the order of one flux 



quantum through each square, in agreement with our theoretical flndings. The tempera- 
ture dependence experimentally observed seems however less drastic than the theoretical 
prediction. 
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VI. GENERIC INTEGRABLE AND CHAOTIC SYSTEMS 



In sections fV] and |V| we have studied in detail specific geometries of conceptual as well 
as experimental relevance. In particular, we have demonstrated the degree of accuracy of 
our semiclassical approach by a careful comparison with exact quantum results. The aim of 
the present section is to take a broader point of view and to give more general semiclassical 
implications concerning the magnetic properties of ballistic quantum dots. We shall first 
consider the weak-field behavior of generic integrable systems, generalizing the results of 
the previous section. We focus on weak fields because only this regime is affected by the 
integrability of the dynamics at zero field. The case of systems which remain integrable at 
arbitrary field strength was discussed in section In the second stage we shall turn to 
chaotic systems (at weak as well as finite fields) and finally finish the section by discussing 
the similarity and differences of the magnetic response for the various cases of classical 
stability. 



A. Generic integrable systems 

We consider the generic magnetic response of two-dimensional integrable systems per- 
turbed by a weak magnetic field breaking the integrability. The Eqs. ( |2.15| ) and ( ^.2U[ ), 
which relate the thermodynamic functions AF^^^ and AF^"^^ to the oscillating part d°^^{E) 
of the density of states, are general relations which apply in particular here. The main 
difficulty is therefore to obtain semiclassical uniform approximations for d°^^{E) interpo- 
lating between the zero field regime, for which the Berry- Tabor Formula |5^ , |53[| (suitable 
for integrable systems) applies, and higher fields (still classically perturbative however), for 
which the periodic orbits which have survived under the perturbation are sufficiently well 
isolated in order to use the Gutzwiller trace formula [^. This problem of computing for 
a generic system the oscillating part of the density of states in the nearly but not exactly 
integrable regime has been addressed by Ozorio de Almeida . We are going to follow 



this approach for the case of a perturbation by a magnetic field. However, for the sake of 
completeness and in order to define their regime of validity, we will give a brief derivation 
of the basic results needed. This is the subject of section |VI A 1| . In section |VI A2| we then 
deduce the grand-canonical and canonical contributions to the susceptibility. 



1. Perturbation theory for magnetic fields 

Let 7i(p, q) be a quantum Hamiltonian which classical analogue can be expressed as 

^(p,q)=^°(p-^A,q) . (6.1) 

7Y*^(p, q) is the Hamiltonian describing the motion in the absence of a magnetic field and 
A is the vector potential generating a uniform magnetic field H. is supposed to be 
integrable which permits to define action-angle coordinates fi,'^2 ^ [0,27r] such that 

at zero field the Hamiltonian 7i*^(/i, I2) depends only on the actions. 
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To compute d°^'^{E) we start from the same basic equations as for the square geometry. 
In the weak-field regime which we are considering, the only recurrent trajectories of the 
sum Eq. ( |5.2| ) which contribute noticeably to the trace Eq. ( |5.1| ) are those which merge 
into periodic orbits of the unperturbed Hamiltonian as H ^ 0. Considering only these 
contributions, which we can label by the topology M of the unperturbed periodic orbits, 
and dropping the Weyl part of the trace Q{E) of the Green function we can write 



g{E) ^ 5] , QUE) = 4^=^ / d^id^2 exp 



2 _ . TT 

n 2 



(6.2) 



Let us now focus on the contribution of the family of closed orbits M. For sufficiently 
low fields we will employ (as in sections |rV| and |V|) that the change in the semiclassical 
Green function by changing H is essentially given by the modification of the phase, Su/'h 
being large in the semiclassical limit. The variation in the determinant Dm can usually 
be neglected. Therefore, in the evaluation of the integral in Eq. ( |6.2|) one should keep the 
(unperturbed) zero'th order approximation for D^^ and evaluate the action up to the first 
order correction. For the action this yields 

'S'M(q,q) = 5^ + 55'M(q,q) (6.3) 

with 

Sl = i p-dq = i I-dip = 27r/M ■ M , (6.4) 

J orbit J orbit 

noting Im the action coordinates of the periodic orbit family M at i7 = 0. The contribution 
(5S'm is expressed in terms of the area enclosed by the unperturbed orbit by means of Eq. (4.7). 



is constant for all members of the family, but 6S generically depends on the trajectory on 
which the point q lies. However, the area enclosed by the orbit and thus SSm. does not change 
when varying q along the orbits. It is therefore convenient to use a coordinate system such 
that one coordinate is constant along the unperturbed trajectory. Writing M = {rui,ru2) 
where Ui and U2 are coprime integers, this is provided explicitly by the standard canonical 
transformation (I, (p) (J, 9) generated by F2(J, (p) = (u2V^i — Uiip2)Ji + V^2<^2 : 

02 = ip2 J2 = h + {Ui/U2)h ' ^ ' ' 



for which 9i is constant along a trajectory on the torus Im- Then 9i specifies the trajectory 
and 6*2 the position on the trajectory. For a square geometry, 6i and 62 are up to a dilatation, 
respectively, the variables Xq and s introduced in section |V|. O2 should be taken in the range 
[0,2^2] (rather than [0,27r]) to ensure that the transformation Eq. ( |6.5| ) constitutes a one 
to one correspondence. 

After substituting q by 6* in the integral of Eq. (|6.2|), 5S depends only on 6'i, but no 
longer on 92- One can moreover show (see appendix 0) the following relation for the zero 
field approximation of the determinant D^: 




^ „ 1/2 ' (6.6) 

2-nrulg"E 
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where I2 = gsih) is the function introduced in section |^ to describe the energy surface E. 
From Eq. ( |6.2|) and (|6.6| ) one gets 



QUE) 



1 1 



1 



2'KrulgE 



1/2 



exp 



TX 



2™2 d6'2 /■^'^ 

J, 



dOi exp 



(6.7) 



The integral over 62 is the period Tm/t of the primitive periodic orbit. In the absence of a 
field the integral over 9i is simply 2tt which gives 



QUE) 



9e 



1/2 



exp I 



h 



TT 



TT 



V^^ 4 



(6.8) 



d^E), the zero field contribution of the orbits of topology M to the oscillating part of 
the density of states, is obtained from Eq. (|6^ ) as d^E) = — (gs/vr)Im^^(i?). Therefore, 
except for the evaluation of the Maslov indices that we have disregarded here, one recovers 
in this way for the integrable limit the Berry- Tabor formula Eq. ([4.3|) of a two-dimensional 
system (as we have used in section |rV|). 

Inspection of Eq. (|6.7|) for weak magnetic fields shows that, upon perturbation, Q^ is 



(6.9) 



just given by the product of the unperturbed result Q^ and a factor 



CUH) = ^ r de, exp 
Zn Jo 



2m 



HA 



$0 



This accounts for the small dephasing between different closed (in configuration space) orbits 
of topology M due to the fact that the resonant torus on which they are living is slightly 
broken by the perturbation. (An orbit of topology M closed in configuration space is then 
generally not periodic, i.e. closed in phase space.) Supposing the unperturbed motion to be 
time reversal invariant, it can be seen moreover that only the real part of C^H) h^is to be 
considered: The function Am.{Oi) is defined for the unperturbed system. Therefore, the time 
reversed of a trajectory labeled by 9i is a periodic orbit of the unperturbed system which en- 
closes an area — ^m(6'i). Its contribution cancels the imaginary part of exp[2i'nHAM.{6i) I ^o]; 
and one can use 



CUH) 



27r Jo 



2tv 



dOi cos 



2tx- 



HA 



(6.10) 



instead of Cm.{H). Since Cm.{H) is real, one obtains from Eq. ( [5.1|) 

d"'\E) = Y: CUH)dl{E) , (6.11) 

Mt^O 

where d^E) is the zero-field contribution given by the Berry- Tabor expression of Eq. ( [4.31) . 
At zero field we obviously have Cm(0) = 1. At sufficiently large field, the integral ( |6l9| ) can 
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be evaluated using stationary phase approximation]^ Cm can be expressed as a sum over 
all extrema of ^m(^i) (^-c of 6S). These are all the periodic orbits which survive under the 
perturbation. It can be seen [|72| that, in this approximation, Eq. ( |6.11| ) yields exactly the 
Gutzwiller trace formula for which the actions, periods and stabilities of the periodic orbits 
are evaluated using classical perturbation theory. Eq. ( |6.11[ ) thus provides an interpolation 
between the Berry- Tabor and Gutzwiller formulae. 

The functions ^m(6'i), and therefore Cm.{H), are system and trajectory dependent. One 
can, however, gain some further understanding of the perturbative regime by following again 
Ozorio de Almeida and writing Aj^{9i) in term of its Fourier series 

oo 

>4m = E'4m sin(n^i-7W) . (6.12) 

n=0 

If is a smooth function of 6i, the coefficients Am^ are usually rapidly decaying functions 
of n. For systems where one can neglect all harmonics higher than the first one, the inte- 
gral Eq. (|6.10|) can be performed, and it is possible to distinguish two types of functions 
Cm.{H), depending on the symmetry properties of the unperturbed family of orbits under 
time reversal. 

Indeed, one may encounter two different situations depending on whether the torus Im 
is time reversal invariant (e.g. square geometry) or has a distinct partner Im* in phase space 
which is its counterpart under time reversal (e.g. circular geometry). In the former case, the 
origin of the angles 6i can be chosen such that ^m(^i) is an antisymmetric function, while 
in the latter case it can be in principle any real function of 6*1.0 

If Jm is time reversal invariant, Am^{—6i) = —Am^{6i) implies that .4m'' = (as well as 
all the phases 7*^"')). In this case 

C{H) ~ Jo(27ri/^W/$^) . (6.13) 

It is interesting to compare the approximation of C{H) given by the above Bessel function 
with the exact integral Eq. ( p.llj ) obtained in section for the shortest family (M = (1,1)) 
of the square geometry. Noting that 9i = eTrxg/a and using Eq. (p.7[), the Fourier coefficients 
All of -^11(^1) are given by 



= <; (n7r)3~^ (6.14) 
I n even . 



16 2 



n odd 



^^To be precise the ratio HA/^q rather than the field must be large. Formally, one has to 
consider not the -ff — > 00 limit, which is incompatible with the classical perturbation scheme, but 
an h {i.e. $0) limiting process, which does not change the classical mechanics. In practice 
this means that the fluxes considered are large on a quantum scale, but still small on the classical 
scale. This is achieved at high enough energies. 

^^Note in the former case Cm = Cm, while in the latter Cm 7^ Cm but ^m + ^m* = Gm^!sa{H) + 
^M*CM*(-ff)- 
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Keeping only the first harmonic of Au{6i) amounts to approximate the function C{!f) of 
Eq. ( |5.11 ) by Jo(32(/?/7r^) which, as seen in Fig. |ll|, is an excellent approximation. 



If the torus Im is not its own time reversal, ^m(^i) is not constrained to be an anti- 
symmetric function, and in particular is usually non zero. Neglecting, as above, all 
harmonics of Am{(^i) except the first gives 

CUH) = cos (|2vr^j Jo (^Svr^j . (6.15) 

If moreover ^m"* ^ .Am'*, then the field oscillation frequency is essentially given by the mean 
area ^m'' enclosed by the orbits of the family while the overall decrease is determined by 
the first harmonic coefficient . The circular billiard can be regarded as a particular case 
where ^m'' is non zero while Am^ as well as all other coefficients vanish. 



2. Magnetic susceptibility for a generic integrable system 

/^From the expression ( |6.11[ ) of the oscillating part of the density of states the contri- 
butions x^^^ ^-iid x*-^"* to the susceptibility are obtained by the application of Eqs. ( |2.15|) 
and ( |2.2(J ), which express AF(i) and AF^^) in terms of d°"'{E,H). Taking twice the field 
derivative according to Eq. ( p. 21 ) and introducing the dimensionless quantities 

{A is the total area of the system) one obtains for the grand canonical contribution to the 
susceptibility 

^ = -24WE ^C(^^) • (6.16) 

If one assumes moreover that there are no degeneracies in the length of the orbits, one has 
for the averaged canonical correction 



x!!l._24.Vi:^^(CXm (6.17) 

Xl m Tm gs 

~ ^i3'M3|,;:(iM)|^''^-^''^- 



The field-dependent component of x^^) for weak fields is given by 

1 



which is always negative. Therefore, for an ensemble of integrable structures the magnetic 
response is always paramagnetic at zero field. We shall come back to this point in the last 
part of this section. 
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B. Generic chaotic systems 



Let us now consider generic chaotic systems, more generally, systems where all the peri- 
odic orbits are sufficiently isolated that the trace of the semiclassical Green function Eq. ( ^.l]) 
can be evaluated within stationary phase approximation. In this case the Gutzwiller trace 
formula provides the appropriate path to calculate the oscillating part of the density of states 
(with or without magnetic field). The Gutzwiller trace formula expresses the oscillating part 
of the density of states as a sum over all [here isolated] periodic orbits t as 



1 Q 

d"'''(E,H) = ydt ; dt(E,H) = ————^ —- cos(-^ - at-) . (6.18) 

t n . ; 7r;irt|det(Mt-/)|V2 ^2' ^ ' 

St is the action along the orbit t, Tt the period of the orbit, Mt the stabihty matrix, at its 
Maslov index, and the number of repetitions of the full trajectory along the primitive 
orbit. All these classical quantities generally depend on energy and magnetic field. If, as 
considered above for the integrable case, one is interested in the magnetic response to weak 
field, one can express dt{E, H) in terms of the characteristics of the orbits at zero field by 
taking into account the field dependence only in the actions. Proceeding in exactly the same 



way as in section [IV A|, i.e. grouping together the contributions of time-reverse symmetrical 



orbits, one obtains the same relation as Eq. (O) |50j62|: 

dt{E,H) = rf°cos 



(6.19) 



d^ is the zero-field contribution of the orbit, obtained from Eq. ( p.l8| ) at H = 0, and At is 
the enclosed area of the unperturbed orbit. In the case of a generic integrable system, the 
zero field regime played a peculiar role: except for the circular and annular geometries which 
remain integrable at all fields, a generic integrable system looses its integrability under the 
effect of a perturbing magneti field. For chaotic geometries on the contrary, the zero field 
behavior is not substantially different from that at finite fields (as far as the stability of 
the dynamics is concerned). Since we are discussing the general semiclassical formalism of 
chaotic systems without referring to specific examples we do not need to restrict ourselves to 
weak fields. Within this generic framework the chaotic geometries have the same conceptual 
simplicity as the systems which remain integrable at arbitrary field studied in section 
Namely Eq. (|6.18|) applies independently of the field, and for derivatives with respect to the 



field one can use 

where At{H) is the area enclosed by the trajectory t at the considered field. Therefore 
the computation of the contribution x^^^ and x*"^^ to the susceptibility follows essentially 
along the same lines as described in section AF^^^ and AF^^^ are given by Eqs. (p. 15 ) 



and ( p.20| ), and to leading order in h the derivatives with respect to the field should be 
applied only to the rapidly varying term. As a consequence, taking twice the derivative 
of the contribution of the orbit t to AF^^^ merely amounts to a multiplication by a factor 
(eA)V(c^)^ yielding 
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^"^-24WE^f4^f (6-21) 



Xl Y n 

where dt is given by Eq. ( |6.18| ). Note that Eq. ( |6.21| ) apphes also to systems which remain 



integrable at all fields provided the Berry- Tabor formula Eq. (O) is used instead of the 
Gutzwiller one. For chaotic as well as for integrable systems, x^^^ can be paramagnetic or 
diamagnetic with equal probability. The response of an ensemble of structures is given by 
AF(2), which can be calculated as a double sum over all pairs of orbits 



Xl V rtr[\det{Mt - I) det{Mt> - /)| V 



2 



A — A' \ ^ fSt — Sf , , vr 



COS [pt + (Jti)- 



A J '""'[-^ I-*--*'. 2, 

(6.22) 



A J \ h V ^2 

Here some remarks are in order. Due to the exponential proliferation of closed orbits in 
chaotic systems off-diagonal terms should be considered at low temperatures since near- 
degeneracies in the actions of long orbits may appear, so that their contributions do not 
average out. However, at sufficiently high temperatures where only short periodic orbits are 
relevant, off-diagonal terms (of orbits not related by time reversal symmetry) are eliminated 
upon averaging. At finite field where time-reversal symmetry is broken (more precisely, when 
no anti-unitary symmetry is preserved) only the terms with t' = t survive the averaging 
process, and (at the order of h considered) x*-^^ vanishes since then A = A'- The origin of 
the weak-field response for an ensemble is a consequence of time-reversal symmetry since 
non-diagonal terms involving an orbit and its time reversal have an action sufficiently close 
to survive the average process but an area of opposite sign. Indeed, assuming (in the weak- 
field regime) an ensemble average such that only diagonal and time reversal related terms 
are not affected, Eq. ( |6.22|) reduces to 



Xl 



At zero field the cosine of the surviving terms in Eq. ( |6.23 ) is one and their prefactors positive. 



This merely refiects that the dephasing of time reversal orbits due to the perturbing magnetic 
field necessarily induces on average a decrease of the amplitude of N°^'^, and therefore by 
means of Eq. ( |2.20c| ) a paramagnetic susceptibility. For extremely large distributions in 
systems size, such as those discussed in section |VF| , even the oscillating patterns of Eq. (|6.23|) 



due to the subsequent rephasing and dephasing of the time reversal orbits contributions 
vanish upon smoothing. In this case, only the paramagnetic response related to the original 
dephasing is observed, and the average susceptibility reaches zero as soon as AnA^H/ $0 is 
of the order of 2tt for all trajectories. 



Magnetization line-shape for chaotic systems 

The expressions we have obtained up to now in this subsection do not require the system 
to be actually chaotic, but only that periodic orbits are isolated. They should therefore be 
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valid also for the contribution of isolated orbits in mixed systems, where the phase space 
contains both regular and chaotic regions. This includes for instance the contributions of 
elliptic, i.e. stable orbits, provided they are not close to any bifurcation and the surrounding 
island of stability is large enough. 

For geometries being actually chaotic it is however possible to proceed further and to 
derive a general expression for the line-shape of the field dependent susceptibility, if the 
temperature is low enough. For temperatures such that the cutoff time Tc of the damping 
factor RT{Tt) is of the order of the period of the fundamental periodic orbits, the average 
susceptibility will be dominated by the shortest orbits, whose characteristics are largely 
system dependent. However, for higher a large number of trajectories will contribute to 



X*^^)^), and a statistical treatment of the sum on the r.h.s. of Eq. ( |6.23|) is possible, yielding 



an universal line-shape for the average susceptibility. For sake of clarity, we discuss here 
only the case of billiard like structures, but the following developments can be generalized 
in a straightforward way to any kind of potentials. 

Two basic ingredients are required here in addition to Eq. ( |6.23|) to obtain the mag- 



netization peak line-shape. The first one is the semiclassical sum rule derived by Hannay 
and Ozorio de Almeida ||7^, which states that in sums like Eq. ( |6.23| ) the two effects of an 
exponential decrease in the prefactors on the one hand and the exponential proliferation of 
orbits on the other hand cancel each other yielding 



6{n-T) 1 
4-|det(M,-/)| r- ^^■'^> 

(Note, that in the above sum the contributions of orbits with number of repetitions > 1 
are neglected.) To be valid, this equation requires that the periodic orbits are uniformly 
distributed in phase space which will only be achieved for sufficiently large r. For billiards 
the periods are given, up to a multiplication by the Fermi speed, by the length of the 
orbits and the periods r in Eq. ( |6.24| ) can be replaced by the lengths L. We call the 
characteristic length for which periodic orbits can be taken as uniformly distributed in phase 
space. Typically, L* is not much larger than the shortest period of the system. 

The second ingredient is the distribution of area enclosed by the trajectories. For chaotic 



systems, this distribution has a generic form [p4| , [74| . Namely the probability P/v(G) for a 



trajectory to enclose an algebraic area G after N bounces on the boundaries of the billiard 
is given by 

Pjv(9) = ^=^ exp ( ^ ) . (6.25) 



This result actually follows from a general argument fTj] which in our case can be stated 
as follows: With a proper choice of the origin, the area swept by the ray vector for a 
given bounce is characterized by a distribution, with zero mean value and a width ctn which 
define the parameter of the distribution Eq. ( |6.25| ). For a strongly chaotic system, successive 
bounces can be taken as independent events, which by means of the central limit theorem 
yield the distribution Eq. ( |6.25[ ). Denoting L the average distance between two successive 
reflections and = o-^lL, this is equivalent to 
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Now -Pl(0) is the distribution of enclosed areas for trajectories of length L, and the above 
equation is valid for L larger than a characteristic value Lg, which again is of the order of 
the shortest closed orbit's length. 

For temperature sufficiently low so that Lc > L^.L^, Eqs (|6.24| ) and (|6.26| ) can be used 
to replace the sum over periodic orbits Eq. ( |6.23|) by the integral 



Performing the Gaussian integral over G, and introducing the dimensionless factor ^ = 
2TrHy/aLLc/^0: one obtains the average susceptibility as 



X'^'n - 96 f ^ ) F(0 (6.28) 



Xl V 

where the function F(^) is defined as 

F(n=/ (l-4eV)exp(-2e^x)c/x ; x = L/L^. (6.29) 

JO V sinh X J 

The quadrature cannot be performed analytically (in a closed expression) for arbitrary ^0, 
but it can easily be calculated numerically. As seen in Fig. 0, F(^) has a maximum at ^ = 
with a half-width ~ 0.252. Expansion of F(0 for small ^ yields F(0 ^ n^/Q - (3C(3) + 
27r'^/15)^^ (where ({x) is the Zeta-function). Denoting A = aiLc/A?, the susceptibility at 
zero field is thus given by 



^ ^-(if = 0) = 167r2A , (6.30) 



Xl 

and the value half-width A$ by 

^ = ^A-'/= . (6.31) 

$0 27r ^ ^ 

The experimental observation of Eq. ( |6.28| ) would be a very stringent confirmation for 
the applicability of the whole semiclassical picture developed here. However, two remarks 
are in order: 

(i) it is experimentally usually rather difficult to make a clear cut distinction between the 
function F(^) we obtained and, say, a Lorentzian shape. Therefore, the temperature depen- 
dence (through Lc) of both the height and, more surprisingly, the width of the magnetization 
peak should be observable rather than the precise functional form of Eq. ( |6.28D . 



^^Using for Rt{L) the asymptotic expression Rt{L) = 2(L/Lc) exp(— L/Lc), valid for L > Lc = 
npvF/TT, yields F(C) = (1 - 5^^) / (l + S,^)^ , but the contribution of the range L < Lc is of the same 
order. 
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The physical picture underlying these results is that at a given temperature, the cutoff 
length Lf. determines the length of the orbits providing the main contribution to the suscep- 
tibility. The smaller the temperature, the larger L^. and the longer the contributing orbits. 
The typical areas enclosed by these orbits thus increase, making them more sensitive to the 
magnetic field and yielding a larger susceptibility at zero field and a smaller width since time 
reversal invariance is more rapidly destroyed. The precise temperature dependence of the 
height and the width (and their relationship, which might be useful when ctl is unknown) 
is given by Eqs. (|6.30|) and ( |6.31|) . 

(ii) It should be borne in mind that Eq. ( ti.28| ) gives only the contribution of the diagonal 
part of x^'^\ but does not take into account the contribution of pairs of orbits which are not 
related by time reversal symmetry. Moreover, the statistical approach used implies that 
fairly long orbits are contributing to the susceptibility, which because of the exponential 
proliferation of such orbits should yield an increasing number of quasidegeneracies in their 
length. Therefore, to smooth out these non-diagonal term, one should a priori require that 
the smoothing is taken on a very large range of {kpo). In practice however, and as will be 
discussed in more detail in |^3[, the smooth disorder characteristic of the GaAs/AlGaAs 
heterostructures for which this kind of experiments are done will actually be responsible 
for the cancelation of the non-diagonal terms without affecting (for small enough disorder) 
the contribution we have calculated^ The effects of non-diagonal terms should therefore be 
noticeably less important in actual systems that it might appear in a clean model. 



C. Integrable vs. chaotic geometries 

The magnetic responses of chaotic and integrable systems have similarities and differences 
with respect to their treatment as well as to the resulting susceptibility. The most remarkable 
similarity is the paramagnetic character of the average susceptibility, while the magnitude 
of this response greatly differs for both types of geometries. Concerning their treatment the 
differences arise form the lack of structural stability of integrable systems under a perturbing 
magnetic field. Indeed, for non generic integrable systems such as the ring or circular billiards 
which remain integrable at all fields, the structure of the obtained equations are, except for 
the use of the Berry- Tabor trace formula instead of the Gutzwiller trace formula, the same 
as those for the chaotic systems. For generic integrable systems however, the breaking of 
invariant tori requires a more careful treatment yielding slightly less transparent, though 
essentially similar expressions. 



Without entering into any details, the reason for this is the following. For smooth disorder, 
one should distinguish between an "elastic mean free path" I, and a transport mean free path It 
which is much larger than /. For small disorder, It can be assumed infinite, but long orbits will 
usually be longer than /. As a consequence, the action of each orbit is going to acquire a random 
phase from sample to sample, which is decorrelated for different orbits, but is the same for time 
reversal symmetric orbits. Thus the diagonal contribution we have calculated will not be affected, 
but non-diagonal terms will be strongly suppressed. 
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1. Paramagnetic character of the average susceptibility 



Because of this formal similarity, the qualitative behavior of the magnetic response is 
also quite the same for generic chaotic and integrable systems. The susceptibility of a 
single structure can be paramagnetic or diamagnetic and changes sign with a periodicity 
in kpa of the order of 2tt. On the other hand, the average susceptibility for an ensemble 
of microstructure is, as expressed by Eqs. (|6.17|) and (|6.23|) , paramagnetic at zero field 
independent of the kind of dynamics considered. Indeed Eq. ( |2.20cj ) states that AF(2) is, up 
to a multiplicative factor, the variance of the [temperature smoothed] number of states for a 
given chemical potential /i. In integrable and chaotic systems the basic mechanism involved 
is that the magnetic field reduces the degree of symmetry of the system, which as a general 
result lowers this variance. Therefore the AF^^) necessarily decreases when the magnetic 
field is applied and the average susceptibility is paramagnetic at zero field. 

There are some differences worth being considered. First, for chaotic systems the only 
symmetry existing at zero field is the time reversal invariance, while for integrable systems 
the breaking of time reversal invariance and the breaking of invariant tori together reduces 
the amplitude of N°'^'^{E). For chaotic systems the paramagnetic character of the ensemble 
susceptibility arises as naturally as the negative sign of the magnetoresistance in coherent 
micro structures. The situation is similar to a random matrix point of view, where the 
ensembles modeling the fluctuations of time reversal invariant systems are known to be less 
rigid (in the sense that the fluctuation of the number of states in any given stretch of energy 
is larger) compared to the case where time reversal invariance is broken. The transition 
from one symmetry class to the other can be understood by the introduction of generalized 
ensembles whose validity can be justified semiclassically [^. It is however important to 
recognize that even for the chaotic case we do not have the standard GOE-GUE transition 



22 1 since (|2.20c|) involves the integration over a large energy interval. We are therefore not 



in the universal, but in the "saturation" regime where {N°^'^{E)y is given by the shortest 
periodic orbits. 

Secondly, for chaotic systems and for temperatures sufficiently low that a large number 
of orbits contribute to the susceptibility, it is possible — similar as in the weak localization 
effect in electric transport p4| — to derive a universal shape of the magnetization peak. 



This is not possible for integrable systems, which do not naturally lend themselves to a 
statistical treatment. 



2. Typical magnitude of the magnetic susceptibility 

Even if there are some analogies between the magnetic response of chaotic and integrable 
systems (especially when the latter remain integrable at finite fields) , the magnitude of the 
susceptibility exhibits significant differences. The contribution of an orbit to the Gutzwiller 
formula for two-dimensional systems is half an order in h smaller than a term in the Berry- 
Tabor formula for the integrable case. More generally, in the case of / degrees of freedom, 
the h dependence of the Berry-Tabor formula is h~^^'^^^^'^ being the same as in the semi clas- 
sical Green function. The Gutzwiller formula is obtained by performing the trace integral of 
the Green function by stationary phase in / — 1 directions, each of which yielding a factor 
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h^^"^. This results in an entire behavior independent of / for a chaotic system. Impor- 
tant consequences therefore arise for the case of two-dimensional billiards of typical size a 
at temperatures such that only the first few shortest orbits are significantly contributing to 
the free energy, and gives rise to a different parametrical kpa characteristic of integrable and 
chaotic systems. The k^a behavior of the density of states and susceptibility for individual 
systems as well as ensemble averages is displayed in Table |. While the magnetic response 
of chaotic systems results from isolated periodic orbits, it is the existence of families of flux 
enclosing orbits in quasi- or partly integrable systems which is reflected in a parametrically 
different dependence of their magnetization and susceptibility on kpa (or \/N in terms of 
the number of electrons). The difference is especially drastic for ensemble averages where 
we expect a kpa independent response x for a chaotic system while the averaged suscepti- 
bility for integrable systems, e.g. the ensemble of square potential wells in the experiment 
discussed in section 0, increases linearly in kpa. Under the conditions of that measurement 
p5| the enhancement should be of the order of 100 compared to an ensemble of chaotic 



quantum dots. We therefore suggested fS^] to use the different parametrical behavior of the 
magnetic response as a tool in order to unambiguously distinguish (experimentally) chaotic 
and integrable dynamics in quantum dots. We stress that this criterion is not based on 
the long time behavior of the chaotic dynamics but on short time properties, namely the 
existence of families of orbits contributing in phase to the trace of the Green function of 
integrable systems. 
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VII. NON PERTURBATIVE FIELDS: BOUNCING BALL AND DE HAAS VAN 

ALPHEN OSCILLATIONS 



Up to now we have essentially focused on mesoscopic effects in the weak magnetic field 
regime where the classical cyclotron radius Tc is large compared to the typical size a of the 
system, i.e. 

^'S>1. (7.1) 



a eH a 

Then, electron trajectories can be considered as straight lines between bounces and the 
dominant effect of the magnetic field enters as a semiclassical phase in terms of the enclosed 
flux. Nevertheless, as shown in Fig. Q in the introduction (for the case of a square) the 
low-field oscillations of x are accurately described by classical perturbation theory in terms 
of the family (11) of unperturbed orbits (left inset in Fig. [l|(b)). They persist up to field 
strengths ^ 10 which is by orders of magnitude larger than the typical flux scale which 
describes the breakdown of first order quantum perturbation theory, i.e., magnetic fluxes 
where the first avoided level crossings appear. Due to condition ( [7. 1| ) the relevant classical 
"small" parameter is H/kp. The semiclassical "weak-field" regime increases with increasing 
Fermi energy. 

In this section we will go beyond this (classically) perturbative regime and discuss mi- 
crostructures under larger fields, where the magnetic response reflects the interplay between 
the scale of the confining energy and the scale of the magnetic field energy %uOc on the quan- 
tum level. Classically, non-perturbative fields affect the motion not only through a change 
of the actions (by means of the enclosed flux), but additionally due to the bending of the 
trajectories. A priori, the semiclassical approach we used for weak magnetic fields applies 
also to this case without any difference: Oscillating components of the single-particle den- 
sity of states can be related to periodic (or nearly periodic) orbits by taking the trace of the 
semiclassical Green function. The magnetic response is then obtained from integration over 
the energy and taking the derivatives with respect to the magnetic field. These operations 
correspond to the multiplication by the inverse of the period of the orbit, by the damping 
factor Rt and by the area enclosed by the orbit. Three field regimes (weak (a ^ rc), inter- 
mediate (a ~ Tc), and high (a > 2rc) fields) can be clearly distinguished as is illustrated in 
Fig. 0(b) for the square geometry. The distinction of the three regimes appears not because 
they deserve a fundamentally different semiclassical treatment, but simply because of some 
salient features of the classical dynamics associated to each of these regimes. 

In the high-field regime, most of the orbits simply follow a cyclotron motion. In that 
case, the system behaves essentially as an infinite system, and one recovers the well known 
de Haas- van Alphen oscillations for x'^^^ We shall moreover see below that within our 
semiclassical approach, the destruction of some of the cyclotronic orbits due to reflections at 
the boundaries can be taken into account, allowing to handle correctly the crossover regime 
where a > 2rc but is not yet negligible with respect to a. 

While the high field (a ^ Tc) classical dynamics is generally (quasi) integrable the 
dynamics in the intermediate field regime is always mixed (in the sense that chaotic and 
regular motion coexists in phase space) except for particular cases of systems with rotational 
symmetry which remain integrable independent of the magnetic field. In contrast to that. 



59 



systems in the small field regime can exhibit any degree of chaoticity in the zero field limit. 
Indeed, there is a large variety of geometries for which the motion of the electrons in the 
absence of a magnetic field is either integrable, or completely chaotic. Therefore, increasing 
the field starting from an integrable (respectively chaotic) configuration at if = 0, the 
intermediate field regime will be characterized by an increase (respectively a decrease) of 
the degree of chaos of the classical dynamics, which will noticeably affect the magnetic 
response of the system. However, if the zero-field configuration already shows a mixed 
dynamics (which is generically the case), the only noticeable difference between the weak 
and intermediate field regime will consist in the complete lost of time reversal symmetry 
and naturally its consequences on x"-^^ as discussed in section [V^. 

In addition, for some particular geometries, namely those for which the boundary con- 
tains some pieces of parallel straight lines, the intermediate field susceptibility will be char- 
acterized by the dominating influence of bouncing-ball orbits, periodic electron motion due 
to reflection between opposite boundaries. Fig. |l](b) depicts a whole scan of the magnetic 
susceptibility of a square from zero flux up to flux = 55 (3rc ~ a). We can see there, 
and we will discuss in detail below, that there are — besides the small-field oscillations due 
to orbits (11) — two well separated regimes of susceptibility oscillations: The intermediate 
field regime (2rc > a) refiects quantized bouncing-ball periodic orbits (second inset) and 
the oscillations in the strong field regime (2rc < a) which, as mentioned above, are related 
to cyclotron orbits (right inset). Although the results to be reported are of quite general 
nature we will discuss them quantitatively for the case of square microstructures. We study 
individual squares and perform our analysis within the grand canonical formalism. 



A. Intermediate fields: Bouncing-ball magnetism 

The full line in Fig. |13|(a) shows the quantum mechanically calculated (see section 0.D) 
grand canonical susceptibility for small and intermediate fiuxes at a Fermi energy corre- 
sponding to ~2100 enclosed electrons in a square at a temperature such that fc^T/A = 8. 
The semiclassical result X(ii) from the family (11) (Eq. ( ^.18|) ) shown as the dashed-dotted 



line (with negative offset) in Fig. |T3|(a) exhibits the onset of deviations from the quantum 
result with respect to phase and amplitude starting at ip ^ 8 ^ 2a) indicating the 
breakdown of the family (11) of straight line orbits. With increasing fiux we enter into a 
regime where the non-integrability of the system manifests itself in a complex structured 
energy level diagram (see Fig. |l|(a)) on the quantum level and in a mixed classical phase 
space of co-existing regular and chaotic motion. However, besides the variety of isolated 
stable and unstable periodic orbits there remains a family of orbits with specular refiections 
only on opposite sides of the square. We will denote these periodic orbits shown in Fig. ^ 
which are known as "bouncing-ball" orbits in billiards without magnetic field by {M^, 0) and 
(0, My) according to the labeling introduced in section 0.A. [M^ and My are the number of 
bounces at the bottom and left side of the square.) These orbits form families which can be 
parameterized, e.g., for the case (M^, 0) in terms of the point of refiection xq at the bottom 
of the square. We thus expect — as in the case of the families {M^, My) in section |V] — 
in the semiclassical limit a parametrical dependence on kpa of the related susceptibilities 
which should strongly dominate the contributions of the co-existing isolated periodic orbits. 
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We present our semiclassical calculation of the susceptibility contribution related to 
bouncing-ball orbits for the primitive periodic orbits, i.e., (Mj.,0) = (1,0) and generalize 
our results at the end to the case of arbitrary repetitions. We proceed as in section |V| for the 
derivation of Xn ■ However, while those calculations were performed in the limit of a small 
magnetic field (assuming i7-independent classical amplitudes and shapes of the trajectories 
(11)) we now have to consider explicitly the field dependence of the classical motion. The 
contribution to the diagonal part of the Green function of a recurring path starting at a 
point q on a bouncing-ball orbit reads 



^io(q, q' = q;^,^^) 



Dio exp 







Hit 





Simple geometry yields for its length, enclosed area, and action 
sm C n 



(7.2) 



k{Lw + 



rc{H) ^ 



(7.3) 



where C, the angle between the tangent to a bouncing-ball trajectory at the point of reflection 
and the normal to the side, is given by (see Fig. 0) 



sinC 



a 

2^ 



(7.4) 



The Maslov index r/io is 4 and will be therefore omitted from now on. 

As in section 0, we will use as conflguration space coordinates the couple q = (xo,s), 
where Xq labels the abscissa of the last intersection of the trajectory with the lower side 
of the square (see Fig. p^ ) and s is the distance along the trajectory. This choice has the 
advantage that Diq{xq, s) is constant, and therefore taking the trace of the Green function 
merely amounts to a multiplication by the size of the integration domain. As discussed in 
more detail in Appendix 0, the semiclassical amplitude Diq is given by 



^io(q,q' = q) = ^ 

\s\ 



dx'n 



dp. 



xo 



XO—XO 



(7.5) 



where (xcp^o) — * i^o^P'xo) Poincare map between two successive reflections on the 

lower side of the bilhard. Noting = {pxo ~ eAx/c)/{hk) {uxq is the projection of the 
unit vector parallel to the initial velocity on the x axis) one obtains from simple geometrical 
considerations 



Xo 
Xn 



Pxo 

Xo 



For the periodic orbits, Xq = xo implies that = a/2rc = sin(, and therefore 



(7.6) 



, 1 hkcosC 

^lo(q,q=q) = ^J^^ 



(7.7) 
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which reduces to Eq. ( p.4| ) in the hmit H = {( = 0). For the contribution of the whole 
family (1,0) we must perform the trace integral Eq. (|5.1| ). The integral over s gives as usual 
a multiplication by the period 



^10 



hk/m 

of the orbit. Moreover, since neither the actions 6*10, nor the amplitude Diq depend on xq, 
the Xo-component of the trace integral simply yields a length factor 



/(i7) = a 1-tan- (7.8) 




(see Fig. |TJ) which describes the magnetic field dependent effective range for the lower reflec- 
tion points of bouncing-ball trajectories (1,0). 1{H) vanishes for magnetic fields correspond- 
ing to 2rc = a. We therefore obtain for the bouncing-ball contribution liio = — (gs/vr)Im^io 
to the density of states 

ME, H) = -j0^,m)L,,D,, sin + ^ • (7-9) 

In order to compute the contribution Xio ^^e (grand canonical) susceptibility we first 
have to calculate by performing the energy integral Eq. (|2.15c| ), and then to take 

twice the derivative with respect to the magnetic field. In a leading % calculation, integrals 
and derivative should again be applied only on the rapidly oscillating part of diQ. Noting 
moreover that Eq. ( ^.7| ) is not restricted to perturbation around H = 0, i.e. that at any field 

dSio _ e 
on c 

we therefore obtain in the same way as we did for Eq. ( |6.21|) 



(1) _ 1 feAio 



2 



rio=-(^] dr,{ii,H)RT{U,) . (7.10) 



Inserting the expressions Eqs. ( |7.3|) , (|7^) and ( [7.71) into Eqs. (|7.9|) and ( [7.10|) , we finally 



have Xio^ explicitly in terms of C, as 

xS_ 3 „^ ^,3 v/^(sinC + cosC-l)(2C-sin(2C))^ 

— - ^2^^-^)' I -4^ X (7.11) 

X sin + ^ RriLio) . 

The entire bouncing-ball susceptibility (xio'' + Xoi)/Xl = '^Xio/Xl according to 
Eq. (|7.11|) is shown in Fig. [l5(a) as the dashed line. At fluxes wp to ip ^ 15 it just ex- 



plains the low frequency shift in the oscillations of the quantum result indicating that the 
overall small field susceptibility is well approximated by Xii + Xio + Xoi- For fluxes between 
ip ^ 15 {vc = 1.2a) up to 9? ~ 37 (the limit where Tc = a/2, i.e., the last bouncing-ball orbits 
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vanish) the magnetic response is entirely governed by bouncing-ball periodic motion and 
the agreement between the semiclassical prediction and the full quantum result is excellent. 

The flux dependence of the actions S'lo (see Eq. ( |7.3| )) is rather complicated. However, 
an expansion for a/vc = 2nip/{kpa) -C 1 yields a quadratic dependence on ip 



^— ~ 2 fcf, a 
n 



1 



1 /27rv9\ 



24 V A; 



(7.12) 



The susceptibility from Eq. (|7. 1 1 ) with S'lo according to Eq. (|7.12 ) is shown as dotted curve 



in Fig. 0(a). It agrees well at moderate flelds and runs out of phase at a flux corresponding 
to a/vc > 1. While the period of the xii small fleld oscillations is nearly constant with 
respect to we flnd a quadratic ip characteristic for the oscillations in the intermediate 
regime which turns into a. 1/ip behavior in the strong fleld regime (see next subsection). 
To show that the agreement between the semiclassical (dashed) curve and the quantum 



result is not an artefact of the particular number of electrons chosen, Fig. |T3|(b) depicts 
semiclassical and quantum bouncing-ball oscillations for fc^T/A = 7 and at a different 
Fermi energy corresponding to ~1400 electrons. With decreasing Fermi energy the upper 
limit Tc = a/2 (or kpa/{27C(p) = 1/2) of the bouncing-ball oscillations is shifted towards 
smaller fluxes {(p ~ 30 in Fig. [T^(b)) and the number of oscillations shrinks. The oscillations 
for (y9 > 30 belong already to the strong fleld regime discussed in the next subsection. 

Up to know we discussed the magnetic response of the family of primitive orbits (1,0) 
and (0,1) which completely describes the intermediate fleld regime at rather high tempera- 
tures corresponding to a temperature cutoff length in the order of the system size. At low 
temperatures we have to include contributions from higher repetitions (r, 0), (0,r) along 
bouncing-ball paths. Lro and Aro have a linear r-dependence, and from the Poincare map 
Eq. ( |7.6|) , one obtains that Dj-o = t~^^^-Dio- Therefore 



(XrO + Xor 



Xl Xl r=i 

^2 



3 „^ v/^(sinC + cosC-l)(2C-sin(2C))- 

YE^'^p^! 7 ^ I'-J-Jj 



SlO TT 

n 4 



r=l 



xY.r-^" sin(r^ + ^) Rrir L,,) . 



Fig. |T^(c) shows the susceptibility at the same Fermi energy as in Fig. |T3|(b) but at a 
signiflcantly lower temperature /c^T/A = 2. The bouncing-ball peaks are much higher 
and new peaks related to long periodic orbits differing from the bouncing-ball ones appear. 
However, the bouncing-ball peak heights and even their shape (which is no longer sinusoidal 
and symmetrical with respect to x = 0) is well reproduced by the analytical sum Eq. ( [7.13| ) 
showing the correct temperature characteristic of the semiclassical theory. 

The kpQ behavior of the bouncing-ball susceptibility at a flxed flux is not as simple as in 
the case of the weak-fleld oscillations (where Xii ~ (^f^)^^^) since the angle C, occurring in 
the prefactor in Eq. ( [7.11] ) depends on kpa and the action is non-linear in kpa. Nevertheless, 
the overall oscillatory behavior is similar as for example in Fig. H(a). However, at a given 
non-zero magnetic fleld the classically relevant parameter Eq. ( [7.1| ) changes with energy. 
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Therefore, by increasing the Fermi energy beginning at the ground state one generally passes 
from the strong field regime (at small energies or high field strengths, see next section) to the 
bouncing-ball regime and will finally reach the regime of oscillations related to the family 
(11). A unique behavior of periodic orbit oscillations is only expected by changing magnetic 
field and Fermi energy simultaneously in order to keep the classical parameter Eq. ( [r.l| ) 
which determines the classical phase space of the microstructure constant. Such a technique 
is known as scaled energy spectroscopy in the context of atomic spectra |7BI . 



Bouncing-ball oscillations are expected to exist in general in microstructures with parts 
of their opposite boundaries being parallel and in spherical symmetrical microstructures as 
the disk discussed in section (In the latter case the oscillations should be even stronger 
than in the square since the effective length 1{H) (Eq. (|7.8|) ) is not reduced with increasing 
magnetic field.) An investigation of rectangular billiards for instance shows a splitting of the 
frequencies of oscillations related to orbits (Mj.,0) and (0, M^y) due to the different lengths 
of the orbits in x and y direction. 



B. Strong field regime 

At large magnetic field strengths or small energy the spectrum of a square potential well 
exhibits the Landau fan corresponding to bulk-like Landau states being almost unaffected 
by the system boundaries, while surface affected states fill the gaps between the Landau 
levels and condensate successively into the Landau channels with increasing magnetic field 
(see, e.g.. Fig. |I](a)). This spectral characteristic corresponds to susceptibility oscillations 
which emerge with increasing amplitude for fluxes corresponding to < a/2, for instance for 
> 40 in Fig. 13(b). They are shown in more detail in Fig. |l^ where the full line depicts the 
numerical quantum result. These susceptibility oscillations exhibit the same period ~ l/i7 
as de Haas-van Alphen bulk oscillations but differ in amplitude, because here the cyclotron 
radius is not negligible compared to the system size. 

For the bulk or in the extreme high field regime Tc <^ a, where quantum mechanically 
the influence of the boundaries of the microstructure on the position of the quantum levels 
can be neglected, an expression for the susceptibility is most easily obtained by Poisson 
summation of the quantum density of states as was briefly sketched in the introduction fol- 
lowing standard textbooks 0]. One then obtains the bulk magnetism as given by Eq. ( |1 .111) . 



It may be interesting to note however that a semiclassical interpretation of this equation 
follows naturally from an analysis similar to the one we followed throughout this paper. In 
this case only one type of primitive periodic orbits exists, the cyclotron orbits with length, 
enclosed area, and action given by 

Lo{H) = 27rr, ; Ao{H) = -nr^^ ; ^ = kL^ + ^HAo = kirr, . (7.14) 

n en 

Moreover, the trajectory passes through a focal point after each half traversal along the 
cyclotron orbit. Therefore, using rjn = 2n for the Maslov indices and omitting the Weyl 
part of G, one obtains from Eq. (|5.2|) a semiclassical expression for the diagonal part of the 
Green function 

G(r,r' = r) = ^ L_ ^(-1)"D„ exp(m7rfcr,) , (7.15) 
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in which the main structure of Eq. (|1 . 1 1|) is aheady apparent. A direct evaluation of the 
amphtude Z)„ in configuration space is however comphcated here by the fact that all trajec- 
tories starting at some point r refocus precisely at r (focal point). Therefore, an expression 
like Eq. (|7.5|) for Dn is divergent and cannot be used. A method to overcome this prob- 
lem by working with a Green function G{x, y]p'^, y') in momentum representation for the x' 
direction instead of G{x,y; x' ,y') is described in appendix 0. It yields (see Eq. ( |E12|) ) 



D„ m 



(7.16) 



Inserting this expression in Eq. ( [7.15|) we obtain the oscillating part of density of states 



<r^\E-H)=Y.dn{E,H) = S!ij^^(-l)"cos(n7rfcr,) , (7.17) 

n n 

from which the de Haas-van Alphen susceptibility Eq. (|1 . 1 1|) is obtained by using 

^(1) ^ 1 (^-^]'Y.d^{^^,H)RT{nL,) . (7.18) 

(with To = Lo/vp) which applies for the same reasons as Eq. ( [7.10|) . 

For an infinite system, this direct semiclassical approach to the susceptibility therefore 
yields the same result as the Poisson summation. For billiard systems, it allows moreover to 
take correctly into account the fact that the trajectories too close to the boundary do not 
follow a cyclotron motion. Indeed, as seen in appendix ^ the contribution of cyclotron orbits 
to the susceptibility Eq. ( |1 . 1 1| ) has to be modified when is not negligible compared to a 
by the introduction of a multiplicative factor s{H). It accounts for the effect that the family 
of periodic cyclotron orbits (not affected by the boundaries) which can be parameterized 
by the positions of the orbit centers is diminished with decreasing field since the minimal 
distance between orbit center and boundary must be at least Vc- One therefore obtains for 
a billiard like quantum dot 



,GC 



oo 



^ = -6s{H) {kpTcf (-1)" RA^nnrc) cos (nvrfc^rj , (7.19) 

XL n=l 



where s{H) is given by Eq. ( [E15| ). In the case of the square we find for the area reduction 
factor 



.iH)^(l-2^) (7.20) 

G being the Heavyside step function. The last cyclotron orbit disappears at a field where 
Tc = a/2, i.e. s{ip) = which happens near ^ 38 in Fig. |l^. There the dashed line 
showing the semiclassical expression ( |7.19| ) is in good agreement with our numerical results 
and reproduces the decrease in the amplitudes of the de Haas-van-Alphen oscillations when 
approaching f{rc = a/2) from the strong field limit. This behavior is specific for quantum 
dots and does not occur in the two-dimensional bulk. Corresponding bulk de Haas-van 
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Alphen oscillations under the same conditions as for the curves in Fig. have (nearly 
constant) amplitudes in the order of xIxl ~ 3000. 

The semiclassical curve which only reflects the contribution from unperturbed cyclotron 
orbits agrees with the numerical curve (representing the complete system) even in spectral 
regions which show a complex variety of levels between the Landau manifolds (see Fig. 
Due to temperature cutoff and since angular momentum is not conserved in the square 
the corresponding edge or whispering gallery orbits are mostly chaotic and do not show up 
in the magnetic response. The strong de Haas-van Alphen like oscillations manifest the 
dominant influence of the family of cyclotron orbits. In related work on the magnetization 
of a (angular momentum conserving) circular disk in the quantum Hall effect regime Sivan 
and Imry observed additional high frequency oscillations related to whispering gallery 
orbits superimposed on the de Haas-van Alphen oscillations. 
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VIII. CONCLUSION 



In this work we have studied orbital magnetism and persistent currents of small meso- 
scopic samples in the ballistic regime. Within a model of non-interacting electrons we have 
provided a comprehensive semiclassical description of these phenomena based on the semi- 
classical trace formalism initiated by Gutzwiller, Balian, and Bloch. We have moreover 
treated in detail a few examples of experimental relevance such as the square, circle and 
ring geometries. 

The global picture that emerges from our study can be summarized as follows. The 
magnetic response is obtained from the variation of the thermodynamic potential (or the 

free energy) under an applied magnetic field and therefore, in a non-interacting model, from 
the knowledge of the single-particle density of states. The semiclassical formalism naturally 
leads to a separate treatment of the smooth (in energy) component of the density of states 
(or its integrated versions) and of its rapidly oscillating part. The former is related to the 
local properties of the energy manifold, while the latter is associated with the dynamical 
properties of the system, more precisely to its periodic (or nearly periodic) orbits. For the 
smooth component we have shown that, despite the leading (Weyl) term in an h expansion 
is independent of the field, higher order terms can be computed and give rise to the standard 
Landau diamagnetism for any confined electron system at arbitrary magnetic fields. In the 
high temperature regime, where the rapidly oscillating component of the density of states 
is suppressed by the rounding of the Fermi surface, the magnetic response reduces to the 
Landau diamagnetism. On the other hand, for the temperatures of experimental relevance 
the contribution coming from the oscillating part of the density of states is much larger than 
the Landau term and dominates the magnetic response. Similarly to the case of diffusive 
systems, the susceptibility of a ballistic sample in contact with a particle reservoir with 
chemical potential can be paramagnetic or diamagnetic (depending on /i) with equal 
probability. The fact that the samples are isolated (with respect to electron transfer) forces 
us to work in the canonical ensemble. Because of the breaking of time reversal invariance 
occurring when the field is turned on, this results, for essentially the same reason as in 
the diffusive regime, in a small paramagnetic asymmetry for the probability distribution of 
the susceptibility of a given sample. For generic integrable systems, this effect is reinforced 
by the breaking of invariant tori, which acts concurrently with the lost of time reversal 
invariance. The asymmetry disappears for a flux A<I> inside the system which is of the 
order of the quantum flux $o at a temperature selecting only the flrst few shortest orbits 
contributions, but may be smaller for lower temperature. Measuring the magnetic response 
of an ensemble of structures with a large dispersion in the size or the number of electrons 
magnifles this asymmetry and yields a total response [per structure] which is paramagnetic 
and much smaller than the typical susceptibility for a flux smaller than A$, and zero 
for larger flux. For ensembles with only microscopic differences between the individual 
structures (i.e. A(fc^a) > 27r, but still Aa/a <^ 1 and AN/N <S 1) further oscillating 
patterns in the average susceptibility should be observed for larger flelds. 

Since the oscillating part of the density of states is semiclassically related to the classical 
periodic orbits, the nature of the classical dynamics quite naturally plays a major role 
in the determination of the amplitude of the magnetic response. Indeed, for a system in 
which continuous families of periodic orbits are present, these orbits contribute in phase 
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to the density of states, yielding much larger fluctuations of the density of states than 
for systems possessing only isolated orbits, and therefore much larger magnetic response. 
Families of periodic orbits are characteristic for integrable systems, while for chaotic systems 
the periodic orbits are usually isolated. This different behavior can therefore be referred to as 
the hallmark for the distinction between integrable and chaotic systems. It should be borne 
in mind however that this difference is due to short-time properties, namely the existence 
or absence of families of orbits, rather than to long-time properties such as exponential 
divergence of orbits. In this respect, some atypical chaotic systems, such as the Sinai billiard 
for instance, may show a magnetic response typical for an integrable system because of the 
existence of marginally stable families of orbits. 

The importance of classical mechanics can be illustrated in the [experimentally relevant] 
case of two-dimensional billiard-like quantum dots in the weak-field regime. If the system 
is chaotic, more precisely if the periodic trajectories are isolated, the typical susceptibility 
scales as {kpa)xL, where kp is the Fermi wave number and a the typical size of the dot. By 
comparison, the typical susceptibility of an integrable system scales with [kpaY^'^XL- This 
characteristic behavior of integrable systems is found in the generic case (like the square) 
where the magnetic field breaks the integrability as well as in the non-generic case (like the 
disk) where the system remains integrable at finite fields. The difference due to the nature 
of the classical mechanics is even stronger for measurements on ensembles of structures 
since one obtains a {kpa)xL dependence for integrable systems and no dependence on {kpo) 
for the chaotic ones. The same parametric dependences are obtained for the persistent 
currents in integrable and chaotic multiply-connected geometries. Therefore, the nature of 
the dynamics yields an order-of-magnitude difference in the magnetic response of integrable 
and chaotic systems, which should be easy to observe experimentally (especially for ensemble 
measurements). Finally, for systems with mixed dynamics, for which the phase-space is 
characterized by the coexistence of regular and chaotic motion, the magnetic response should 
be dominated by the nearly integrable regions of phase-space. This gives rise to a {kpaY^'^XL 
dependence for the typical susceptibility as long as some families of periodic orbits remain 
sufficiently unperturbed. The precise calculation of the prefactor may however present some 
complications that we have not considered here (the general semiclassical treatment of mixed 
systems remains an open problem) and should depend on the fraction of phase-space being 
integrable. 

The semiclassical approach we are using not only allows a global understanding of the 
magnetic response of ballistic devices, but also provides precise predictions when specific 
systems are considered. The detailed comparison between exact quantum calculations and 
semiclassical results for the square geometry demonstrates indeed that the semiclassical pre- 
dictions are extremely accurate. This has been shown in section for weak fields, such that 
the trajectories are essentially unaffected by the magnetic field, and also in section |VII| for 
fields large enough to yield a cyclotron radius of the order of the typical size of the structure 
(where the bending of the classical trajectories has to be taken into account). For intermedi- 
ate fields we have identified a new regime where the magnetic susceptibility is dominated by 
bouncing-ball trajectories that alternate between opposite sides of the structure (enclosing 
flux due to their bending). For high fields the electrons move on cyclotron orbits and we 
have recovered the de Haas - van Alphen oscillations (with finite-size corrections that we 
calculated semi classically) . 
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In order to understand the success of the semiclassical approach, it should be kept in 
mind that the lack of translational invariance characteristic for the ballistic regime, where 
the shape of the device plays an important role, complicates the application of other ap- 
proximation schemes as, e.g., diagrammatic expansions. Therefore, except for very specific 
cases where exact quantum calculations are possible, and unless one is satisfied by direct 
numerical calculations, some semiclassical ideas have to be implemented to deal with such 
problems. Moreover, from a more practical point of view, the semiclassical trace formalism 
we have used appears perfectly adapted to deal with thermodynamic quantities such as the 
Grand Potential or its first and second derivatives iV(/i) and D{fi). Indeed, the beauty 
of this approach is that the oscillating part of the density of states is directly expressed 
in terms of Fourier-like components, each of which is associated with a periodic (or nearly 
periodic) orbit. The thermodynamic properties are obtained from their purely quantal (or 
zero temperature) analogues u, n and d by temperature smoothing, which merely amounts 
to multiply each oscillating component by a temperature-dependent damping factor. For 
all fields (high, intermediate, or weak), this factor depends only on the ratio of the period 
T of the corresponding orbit and the temperature-dependent cutoff time Tc = jSh/n and 
suppresses exponentially the contribution of orbits with period longer than Tc. As a conse- 
quence, not only the effect of temperature is taken into account in an intuitive transparent 
way, but in addition only the shortest periodic orbits have to be considered in the semi- 
classical expansion. All the problems concerning the convergence of trace formulae and the 
validity of semiclassical propagation of the wave function for very long times are of no im- 
portance here. One therefore avoids most of the problems which plague the field of quantum 
chaos when semiclassical trace formulae are used to resolve the spectrum on a mean-spacing 
scale. Mesoscopic physics is usually concerned with the properties of the spectrum on an 
energy-scale large compared to the mean-spacing. In the spirit of the work of Balian and 



Bloch this is the situation for which the semiclassical trace formalism is especially 
appropriate. 

Having stressed the success of the semiclassical approach in dealing with our model of 
non-interacting electrons evolving in a clean medium, it is worthwhile to consider in more 
detail how the above picture should be modified when going closer to the real world, and 
incorporating the effects of residual disorder, electron-electron or electron-phonon interac- 
tions. As stressed in the introduction, the first of these points is relatively harmless because 
of finite temperature smoothing. The restriction to short periodic orbits actually justifies an 
approach to the ballistic regime using a model for clean systems since long diffusive trajec- 
tories do not contribute to the finite-temperature susceptibility. Indeed, careful numerical 
and semiclassical studies of the effect of small residual disorder show that, except for a 
possible reduction of the magnetic response, the above description of the orbital magnetism 
of ballistic systems remains essentially unaltered. In particular, the mechanism proposed 



by Gefen et al. |Q is not borne out by the numerical simulations at the temperatures of 
experimental relevance. For smooth disorder, such as presumably prevails in the systems of 
Refs. and |^, the magnetic response is decreased by the dephasing of nearby trajec- 
tories in a way that depends on its strength and the ratio between the correlation length 
and the size of the structure but diffusive trajectories can be seen to be absolutely 



irrelevant if the elastic mean free path is larger than the size of the structure. The precise 
knowledge of this reduction is however needed in order to make a decisive comparison with 
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the experimental results of Ref. ||2^ . 

At the low temperatures of the experiments the inelastic mean free path of the electrons is 
much larger than the system size since electron-phonon interactions are suppressed. On the 
other hand, the effect of electron-electron interactions on the magnetic response is a much 
more controversial point. In particular, it has been invoked to be the necessary mechanism 
to obtain the measured values for the problem of persistent currents in disorder metals. 
In a first approximation to the experimental conditions that we investigated in this work we 
would infer that electron-electron interactions are not crucial since the screening length is 
much smaller than the size of the samples and since the 2-d renormalization of the effective 
mass at these electron densities is only about 11% |7^. Clearly the two previous criteria will 
not be satisfied in smaller structures, and the possibility that electron-electron interactions 
express themselves through a mechanism for which these estimates are not relevant remains 
open even in the experimental realizations we consider. 

Contrarily to the effect of disorder, which can be implemented within a semiclassical 
framework without essential difficulties, a semiclassical treatment of the electron-electron 
interaction still remains an open problem. However, the genuine effects that we have found 
within our semiclassical approach for the clean model of non-interacting electrons should 
prevail in more sophisticated theories. We think that the rich variety of possible experimental 
configurations for ballistic devices (the shape and the size can nowadays be chosen at will) 
provides an ideal testing ground for these more complete approaches. We hope that the work 
presented here will stimulate experimental and theoretical activity addressing the magnetic 
response of ballistic microstructures. 



ACKNOWLEDGMENTS 

We acknowledge helpful discussions with H. Baranger, O. Bohigas, Y. Gefen, 
M. Gutzwiller, L. Levy, F. von Oppen, N. Pavloff, B. Shapiro, and H. Weidenmiiller. We 
are particularly indebted to H. Baranger for continuous support and a careful reading of the 
manuscript, and to O. Bohigas for forcing us not to stop until getting to the bones of the 
problem. We thank B. Mehlig for communicating us Ref. KR and RAJ acknowledge 
support from the "Cooperation CNRS/DFG" (EB/EUR-94/41). KR thanks the A. von 
Humboldt foundation for financial support. The Division de Physique Theorique is "Unite 
de recherche des Universites Paris 11 et Paris 6 associee au C.N.R.S.". 



70 



APPENDIX A: CONVOLUTION OF A RAPIDLY OSCILLATING FUNCTION 
WITH THE DERIVATIVE OF THE FERMI FUNCTION 



When considering thermodynamic quantities related to the oscillating part of the density 
of states at finite temperature T, one has to evaluate integrals of the form 



J(T) 



dE A{E) exp 



f{E - /i) 



(Al) 



where f'{E — fi) is the derivative of the Fermi function 

f{E - fi) = — ]r— , 

1 + exp p{E — n) 

and (3 = 1/kBT. The rapidely oscillating function A{E) exp ■^S{E) usually originates from 
the contribution of a classical orbit (or a family of orbits) to the oscillating part of the 
density of states. In this case S{E) is the action integral along the orbit, and its derivative 
dS/dE = t{E) is the period of the orbit. 

At zero temperature /' = —6{E — fi) giving for Iq = /(T = 0) 



Jo = -A{n) exp 



(A2) 



In this appendix we show that, to leading order in h and in /3 ^ (but without making any 
assumption concerning their relative value) , the integral of Eq. ([Al|) is given by 



/(T) = loRrir) 



with the temperature dependence 



RAr) 



sinh(r (/i)/rc 



TT 



(A3) 



(A4) 



For systems without potential, i.e. free particles confined in a box (billiards), the period of 
the trajectory is related to its length L by r(/x) = L/vp, where Vp = Tikp/m is the Fermi 
velocity. Rt can then be written as 



Rt{L) 



L/L, 



hvp(3 



sinh(L/Lc 



TT 



(A5) 



In the case of unconfined free particles, the formulae (|A4|) and ( |A5D are equivalent to the 
usual form of the temperature dependence of the de Haas-van Alphen effect Eq ( p..lO| ) given 
in the introduction. Below we present a slight variation of a standard calculation (see e.g. 0) 
of the temperature dependence of the de Haas- van Alphen effect, which generalizes it to 
any type of dynamics, once we caste it in the form of Eq. (|A4|) . 

Performing the integral ( |Al| ) along the contour shown in Fig. |16] and noting that the 
singularities of the derivative of the Fermi function are double poles located dX E^ = ix + 
i{2k + l)vr//9 {k = 0, ±1, ±2, . . .) with a coefficient one finds the following relation 
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dEA{E)exp 



/•oo+j27r//3 

f\E-^) - / dEA(E)exp 

Jo+i2TT/l3 



2m iT{E 



(3 h 



-A{Ei)exp 



(A6) 



At low temperatures, the function f'{E — fi) is essentially zero everywhere in the complex 
plane, except for a narrow band of width j3~^ near the line Iie{E) = fi, therefore the vertical 
portions of the contour {E = and E ^ fi) give negligible contributions. Noticing that 
f'{E — fi) has a periodicity of 2m / (3 we can ignore the complex part of E in the factor /' 
of the second integral. Finally, since Im(i?i) = TT/f3 and Im(£') = 271/(3 along the upper 
portion, we can evaluate the prefactors at fi and expand the actions (which are multiplied 
by 1/fi) as S{E) = S{fi) + T{fi){E — fi) in leading order in and h, obtaining 



/(T) 1 - exp 



27rr(/i) 
(3n 



27rr(/i) 
(3h 



A{fi) exp 



That is, 



J(T) (^1-exp 



2r(/x) 



r 2r(^) 

Jo exp 



r(/i) 



from which one readily obtains the result of Eq. 



Further comments 

We would like to use the above calculation to motivate some choices made in section |T| 
which might have appeared rather arbitrary. Concerning for instance the grand potential 
two equivalent expressions have been introduced: The usual Eq. ( p..3| ) and Eq. ( p.6c|) 
which is obtained from integration by parts. On the other hand, we have used only 

roc 

= - / dE uj°''iE) f{E - fi) (AT) 
Jo 

as the "operational" definition of the oscillating part of fi, and one might wonder 



whether an integral analogue to the one of Eq. ( |1.3|) like 



- j (r'\E)f'^\E - fi)dE , (A8) 

(where f^^^\E — fi) = In (1 + exp [(3{fi — E)]) / (3 is the primitive of the Fermi function) could 
not be used as well. This is not the case for the two following reasons: 

i) First, the oscillating functions d°^'^ or uj°^^ are usually obtained in a semiclassical ap- 
proach and are therefore valid only for large energies. If the chemical potential fi is 
in the semiclassical regime and (3~^ <^ /i, which is always the case for the problems 
we consider, only the neighborhood of fi in which uj"^'^ can be used safely, contributes 
significantly to the integrals of Eq. ([A^). On the contrary, the integral ( [A^ ) involves 
energies close to zero. Therefore there is no reason that dP^'^ is accurate, being quite 
often a diverging function. 
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ii) In addition, even if one has at hand an equation as ( |1.6| ) which is a non-semiclassically 
exact expression, the integrals of Eq. (^) and Eq. (|A8|) are, strictly speaking, not 
equivalent. The latter form contains some boundary terms not present in the former, 
which obviously have to be removed from Q'^^'^ since they do not average to zero under 
a local smoothing. 

In a semiclassical treatment the derivative /' of the Fermi function is superior to any of its 
integrated versions since it is significant only at energies where semiclassical approximations 
can be used safely. 



APPENDIX B: SEMICLASSICAL EXPANSION OF THE MEAN DENSITY OF 

STATES 

In this appendix we calculate the first two terms in an h expansion of the smooth part of 
the density of states d{E). We follow the standard approach introduced by Wigner using 
the notion of the Wigner transform of an operator. The Wigner transform of a quantum 
operator O is defined by 

[(^]w(q,p)= /dxe-P-/^(q+||a|q-|). (Bl) 

Among different properties of the Wigner transform we will essentially make use of the 
following two: First, the trace of an operator is related to the integral over phase space of 
its transform by means of 

Tr(a) = J dqdp [0]w(q, p) • (B2) 

(We stress that this is an exact, not semiclassical, relation.) Secondly, for any operator 
function of the position and momentum quantum operators J^(q, p) (with some specified 
ordering), the semiclassical leading order approximation to its Wigner transform is just the 
related classical function, that is 

[^(q,p)]w = ^(q,p) + 0(;i). (B3) 

When !F depends only on q or p, the relation between the Wigner transform and the classical 
function is exact (no corrective terms in h), as can be directly checked from Eq. ([BT|). 



We will follow closely the presentation of Ref. [^] to which the reader is referred to for 
further details. The first step in the calculation of d{E) is to consider the Laplace transform 
of the level density, Eq. (p.2|), which due to the property ([B2|) can be written as 



Z{X) = I dqdp [e-^«]w(q, p) (B4) 



(gs = 2 is the spin degeneracy factor). Using Eq. (|B3|) , the leading order [Weyl] term in h, 



Zw, is obtained by replacing [e '^''^]>v(q) p) by e -^^(I'P), where 
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is the classical Hamiltonian. At this level of approximation Z{X) is given by 

(p-fA)2 



^w(A) 



dqdp exp —A 



2m 



V(q) 



(B6) 



Since this term is field independent (see the change of variable (|3.5|) in the text) we need 
to go to the next order in A in order to obtain non-vanishing contributions to the magnetic 
susceptibility. Therefore we consider the asymptotic semiclassical expansion 



[e-^^]w(q,p) = e 



-A[W]w(q,p) 



^^Y Cn{q,p, A) 



E 

n=0 



(2n)! 



(B7) 



We have already seen that Cq = 1. The following coefficients C„ can be obtained recursively 
by grouping terms according to their power in h. In particular, the first coefficient is given 
by MM 



dCi 



A e 



-A[H]v 



(B8) 



where A= E ^li^Pi ~ dpidqi is the Moyal bracket. 



i=l 



Inserting the classical Hamiltonian ( |B5| ) at the place of [H]yv we obtain Ci as the sum 
of three terms 



where 



(^1°) = - 



m 



c, = (c?) + (Q) + (CD 



A^ VV(q) - — (mVV(q) ■ VV(q) + ((p - f A) ■ V)^ V{q] 



(B9) 



is, up to the change of variable Eq. (|3.5|) in the integration over the phase space, the first 
order correction without magnetic field Ci{h=o) given in |]65|. {C^) is given by terms being 
antisymmetric in p that vanish when taking the trace over phase space. Finally, 



2 d 

(CD = -T. 

i,k=l 

2A3 



A^ 



[{dqiAk){dqiAk) - {dqiAk){dqkAi)] + 



+ 



3m^ 



{dqiAk){dqkAi){pi - lA,f - {dqiAk){dqiA^ 



{Pk - ^^Akf + {pi - lAi 



(BIO) 



The first correction to the Laplace transform of the density of states is given by 



Zi(A) 



dqdp 



A^ 



6 



-AW 



(Bll) 



if^B = (e^)/ (2mc) is the Bohr magneton). Eqs. (|B11D and ( p.7|) are obtained from Eq. ( piO|) 
by using the identity 
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E 

jk 



dAj dAk 
dqk dqj 



(B12) 



and a few transformations that leave the integral over p unchanged, namely : (i) the change 
of variables Eq. (|3.5| ) (allowing the substitution of Ti. by Ti" = p^/2m+\^(q)) (ii) the elimina- 
tion of all terms antisymmetric in p, (iii) the replacement of all terms of the form pfe~^^ Z^*" 
by (m/A)e^'^P Z^"^. Note finally that the field appears only in the term — which 
is independent of the confining potential V^(q). This is at the root of the very general 
applicability of the Landau result. 



APPENDIX C: CALCULATION OF qe FOR A RING BILLIARD 

In this appendix we derive the explicit form I2 = gE{h) of the energy surface E in 
action space for a ring geometry. The calculation reduces to the evaluation of the integral 
of Eq. ( |4.1| ) along two independent paths on the invariant torus. The only subtlety arises 
from the difficulty of visualizing the integration paths in our four-dimensional phase-space 
where the tori are discontinuous due to the presence of hard walls. We closely follow the 
procedure used by Keller and Rubinow ||6^ for the circular billiard, and we refer to this 
work for further details. 

In a circular ring (with outer and inner radii a and b) we can distinguish two types of 
periodic trajectories: those which do not touch the inner disk (type-I, Fig. p!7|a) and those 
which do hit it (type-II, Fig. p!7|b). Type-I trajectories have their caustics outside the inner 
disk and therefore they are unaffected by their presence. They have an angular momentum 
pc, with b < c < a {p = \/2mE). Taking as the integration path Ci the concentric circle of 
radius i?, we have pdq = {pc/R)dq and then 

1 f c 

h = — f ~5 P^(l = Pc ■ (CI) 

ZTT JCi R 

The action variable Ji is just the angular momentum. The straight part of the path C2 of 
Fig. p!7| a is chosen along a classical trajectory, where p is constant and coUinear with dq. For 
the part along the outer circle pdq = — (pc/a)dg. Combining both contributions we have 



t{I) _ P 
-'2 — ~ 



TT 



2 2 

a — c 



^ — c arccos \ -] \ ■ (C2) 



Elimination c between (CI) and (|C^ ) leads to Eq. ( 4.12| ) of the text, valid for the description 



in action space of the energy surface of the circular billiard [^] and the energy surface 
associated with type-I trajectories in the ring billiard {pb < Ii < pa). We have chosen 
the integration paths for type-I trajectories different from those of Ref. because slight 
modifications of them are applicable for type-II trajectories. 

Type-II trajectories have their caustics in the interior of the inner disk, that is, they have 
an angular momentum pc, with c < b. Integration along the path Ci of Fig. p!7| b leads to the 
identification of Ji with the angular momentum pc (similarly to Eq. (PT|)). By choosing the 
path C2 as shown on Fig. p!7|b, the action integral along this path is simply the difference 
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ril), 



l2^\b) (where both terms are given by Eq. (C2), except that for the second a should 



be replaced by b). This yields 



r(//) _ P 



2 2 

a — c 



1/2 



1/2 



— C 



arccos 



— arccos 



(C3) 



Eliminating c between (PI]) and (CSl) leads to Eq. ([4.21|) of the text. 



APPENDIX D: CALCULATION OF THE DETERMINANT Dm AT ZERO FIELD 
FOR A GENERIC INTEGRABLE SYSTEM 



In the semiclassical approximation of the Green function (Eq. 
associated with a classical trajectory t is given by [p^ 



]2[)) the amplitude Dt 



dqdq' dqdE 



dEdq' dEdE 



1/2 



,. .,,1/2 

mq2\ 



dqidq'i 



1/2 



(Dl) 



Note the second equality holds not only when q2 is taken along the orbit and qi in the 
perpendicular direction, as supposed by Gutzwiller in its original derivation but also, 
as shown by Littlejohn, in any coordinate system 



(see Sec. III.C in Ref. [^] and Sec. Ill 
in [81 1). Although a priori qi and q2 play a similar role, their non-symmetrical appearance 
on the right hand side of Eq. ( pi|) (gi and q2 can be exchanged without affecting the value 
of -Dm) is due to the fact that one coordinate (here gi) is chosen as a Poincare surface of 



section, and the dependence of the other coordinate (here ^2) just expresses the conservation 
of energy. 

Turning now to the particular problem we are concerned with, i.e. an integrable system 
at zero field and the diagonal part of the Green function, the above Eq. ( [PT| ) applies to Dm 
(except for a change t — > M in the label of the orbits). Moreover, the measure DM.dqidq2 in 
Eq. ( |6.2|) is invariant under the transformation (51,52) (^1)^2) at zero magnetic field (see 
Ref. |]72| for a more detailed discussion of this point). In other words, noting in Eq. ( pi| ) 
Dm{^) the determinant in the original q coordinate and -Dm(^) the determinant defined in 
the same way but in the system of coordinates given by Eq. ( |6.5D , one has DM.{^)dqidq2 = 
D-M.{0)d6id62- Therefore 



I^M(q) 




6262 



, 1/2 



d9id9[ 



1/2 



dJi 



-1/2 



(D2) 



where the derivatives have to be taken at E, 62, and 



02+211 M2 constant. To compute the 
r.h.s. of Eq. (P2| ) one just needs the expression of the Poincare mapping (6'i, Ji) — > {9[, J[) 
between the two {62 = const.) Poincare surfaces of section. Since the motion is integrable, 
J[ = Ji, and from Eq. (|6.5[) we obtain 



^;(Ji, 61) =61 + 2nrul{a{Ji) - M1/M2) 



(D3) 
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where a(Ji) is the winding number of the torus labeled by Ji. Thus 



2'Krut 



da 
dJi 



(D4) 



We recall that the function qe introduced in section |^ is defined by the implicit relation 
H{Ii, l2 = gE{h)) = E, which yields after differentiation dgE/dli = —a. Therefore 



da 



da 



dJi ^ dli 



d?gE 



from which one finally obtains 



9. 



2nrulg"E 



1/2 



(D5) 



(D6) 



Z^M is inversely proportional to the square root of the curvature of the line H{Ii, I2) = E 
const, and independent of 6. 



APPENDIX E: DIAGONAL PART OF THE GREEN FUNCTION FOR A FREE 
ELECTRON IN A CONSTANT MAGNETIC FIELD 

In this section, we calculate semiclassically the diagonal part of the Green function G(r, r) 
for a free electron moving in a plane in a perpendicular magnetic field. The resulting classical 
cyclotron motion is extremely simple, but yields slight complications in the semiclassical 
evaluation of the diagonal part of the Green function because all trajectories starting at 
some point r are refocused precisely at r. The calculation of the prefactors deserves special 
attention but can be done using a slight variation of the standard techniques and yields 
for unconstrained systems the usual result Eq. ( |1 . 1 1| ) . In addition to provide an alternative 
(semiclassical) derivation of the de Haas-van Alp hen effect, our procedure allows to compute 
correctly the contribution of the cyclotron orbits for billiard systems, i.e. it takes into account 
the corrections due to the boundaries which appear to be necessary if the cyclotron radius 
Tc is not small compared to the typical dimension a of the system. 



a. Computation of the prefactor of a Green function near a focal point 

It is an old problem to obtain a correct semiclassical solution of wave equations valid also 
near turning points, focal points, caustics, etc., where the usual expressions are diverging. 



A general solution for this problem can for instance be found in the book of Maslov |^ 
In this subsection, we will give the explicit form of this general theory when applied to the 
calculation of a two-dimensional Green function and consider in the next subsection the 
particular problem of a free electron in a constant magnetic field. To avoid confusion we will 
slightly modify our usual notations, writing G(r|r') instead of G(r,r'). In addition, we will 
make more explicit what are the initial and final points by using r* = {x\y^) for the initial 
(source) point and, and r-^ = {x^ ,y^) for the final (observation) point. As already stated. 
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the semiclassical evaluation of the Green function G{r^\r-^) yields a sum over all classical 
trajectories t joining to r-^ at energy E 



G'(r>0=EG'i(r>^) 



(El) 



For a trajectory t starting at such that r-'^ is not a focal point, one can use (cf. Appendix 0, 
and in particular the discussion concerning the non-symmetric role of x and y) 



1 1 



iJi \/2mh 



\yy 



/1 1/2 



Dt exp 



dx^dxf 



—bt — m — 
h '2. 

1/2 



(E2) 
(E3) 



However, the above expression is not valid near focal points where [locally, and at fixed 
y^] x^ becomes independent of p].. The use of the action integral S{r\r^) supposes that r* 
and r-'^ can be taken as independent variables, and d^S/dx^dx-^ is a priori not meaningful 



— oo one sees 



since x^ is entirely determined by x\ Writing d'^S /dx'^dx^ = —{dx^/dp^^) 
moreover that Dt is, as mentioned above, in fact diverging. 

To overcome this difficulty Maslov proposed a procedure to compute Gt(r*|r-^) using 
a momentum (or mixed position/momentum) representation, by defining (omitting for a 
moment the source point r*) 



G,{xf ,yf) = :F-^^jGU.y^)] 



where T } .is the inverse Fourier transform 

ni. — >xJ 



/[•] 



dpll 



■J exp ( -x^pI 



(E4) 



(E5) 



performing quantum mechanically the change from the mixed representation {p^, y) to the 
position representation (x, y) . 

Eq. (0) is just the definition of Gt which remains to be evaluated semiclassically. The 
general theory presented in Ref. (Section 5.1) can be however applied to our problem, 
giving 



1 



where 



y/2iTrh 
S = S — p{.x^ , 



Dt exp 



Dt = Dt 



dx^ 



dpi 



1/2 



rjt if dpl/dx^ > 
7]t + I if dpi I dx-f < 



(E6) 

(E7) 

(E8) 
(E9) 
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Without entering into a derivation of this semi classical formula for Gt, it can be checked 
that starting from Eq. (^) the evaluation of the inverse Fourier transform Eq. (|E4|) using 
stationary phase approximation readily yields Eq. (|E2|) . Far from any focal point both 
expressions are equivalent at the semiclassical level. Near a focal point however, Eq. ( [E6|) 
still provides an accurate approximation for Gt because the Lagrangian manifold, on which 
the Green function is constructed semiclassically, has a non singular projection onto the 
plane {px,y)- Therefore, contrarily to Eq. ( |E2| ) which is diverging, Eq. ( [E^ ) is still a valid 
semiclassical approximation for Gt, provided the inverse Fourier transformation is evaluated 
exactly (or using uniform techniques going beyond stationary point approximation [^]). 
/^From Eqs. (|E3|) and (pH]) one has 

li/r I 

This explains why the Legendre transform St of St has to be understood as a function of 
and pI- In practice, this means that, to compute St from Eq. (^Tf), the action integral St 
has to be calculated for a trajectory starting at position and ending with a momentum 
pI, and that in the additional term x^pl, x^ has to be interpreted as x^ Finally, 
note that for a Hamiltonian, which can be decomposed into a kinetic energy plus potential 
part (including the case where a magnetic field is present), dp^/dx^ is always negative just 
in front of a focal point and always positive directly after the focal point. Therefore 




rjt right after a focal point 
rjt + 1 just before a focal point 



Since precisely at focal points rjt is incremented by one unit (for a kinetic energy plus 
potential Hamiltonian), this implies that, when crossing a focal point, fjt remains constant, 
keeping the value which rjt acquires after the focal point. This latter has to be taken into 
account for the computation of Maslov indices of a trajectory at a focal point. 



d^S 



dx^dpl 



(ElO) 



b. Application to the cyclotron motion 

Turning now to the specific problem we are concerned with, i.e. cyclotron motion and 
diagonal elements of the Green function, we need to calculate S{x\y\pl,y'l'), where we can 
however restrict ourselves to y^l" = y^ since the partial derivatives are taken only in the x 
direction. Eq. (^) states that, omitting the ?/'s, Sn{x\pl) = Sn{x\p{.) —x-l" {x'' ,pl)pl, where 
Sn{x\pl) is the action integral along a trajectory starting at the abscissa x* and arriving 
with a momentum pi. But here the Poincare map (x*,p^) (x'^,p^) is just the identity, 
and therefore 

Sn{x\pi) = nSo-x'pi , 



dx^dpi 
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where 5*0 is given by Eq. ( |7. 14|) . Noting moreover that ?/* = ij^ for all trajectories, and that 
they pass through two focal points at each turn, one has from Eq. (|E4|) (omitting the Weyl 
part of G) 



G(r*|r^ 



E 



1 



exp 



h 



dpi 



(Ell) 



At fixed position, dp^. = mdx^ (the vector potential eliminates). Therefore the remaining 
integral in Eq. ( [Ellj ) is just, up to a multiplication by the mass m of the electron, an integral 
over the angle 6 specifying the direction of the trajectory at r. For unbounded motion, it 
simply gives a factor 27rm, yielding the expected result 



TTi 

^(-l)"-2exp(mSoA-m/2). 



(E12) 



In billiard systems the contribution to G of the cyclotron orbits is the same as for the un- 
bounded motion, except that for points r close to the boundary, Eq. (|E12|) has to be reduced 
by a multiplicative factor 6'eff/(27r), where 9es is the angular measure of the trajectories not 
affected by the boundary. The contribution to the density of states of the cyclotron orbits 
is thus 



d'^'^E, H) = s{H) y (-1)" cosimrkr. 



(E13) 



The multiplicative factor s{H) is given by 



s{H) 



drd6({r, 



(E14) 



The function ({r, 9) is defined such that C = 1 if the trajectory started at r with initial 
velocity along does not hit the boundary, and C = otherwise. Substituting in the 
integral above the variables (r, 9) by (f , 9), where f specifies the center of the cyclotron orbit 
and 9 the position on this orbit (the Jacobian of the transformation is equal to one) and 
performing the integral over 9 since then C, depends only on f , one obtains 



s{H) 



1 

I 



rfrC(f) 



(E15) 



which yields Eq. (|7.2CI|) for the square geometry. 

As a final comment, we note that the approach described here for a two-dimensional 
electron gas can be generalized in a straightforward manner to three dimensional systems, 
including cases with non-spherical Fermi surfaces. 
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FIGURES 



FIG. 1. (a) Evolution of the first 200 energy levels (of one symmetry class (see section VD| )) of 



a square billiard in a uniform perpendicular magnetic field H as a. function of the normalized flux 
if = Ha? /(^Q ($0 = hc/e). The energies are scaled such that the zero field limit gives E = n'^ ~^^'y- 
At high fields the levels converge successively to the Landau levels while the non-integrable inter- 
mediate field regime exhibits a complex spectral structure, (b) Full line: Numerically calculated 
susceptibility of the square at finite temperature at an energy corresponding to ~1100 occupied 
independent particle states. The susceptibility, being strongly enhanced with respect to the bulk 
value xl^ exhibits pronounced oscillations which are accurately reproduced by analytical semiclas- 
sical expressions (dashed line) based on families of flux-enclosing electron orbits (shown in the 
upper insets for the different magnetic field regimes). 

FIG. 2. Schematic illustration of the separation of the density of states D{^) (solid line) into 
a smooth part D (dashed line) and an oscillating component. The total number of electrons N is 
indicated by the shaded area, and equal to the product of D and p.. 

FIG. 3. Magnetic susceptibility at zero field for a circular billiard of radius a as a function of 
kpa (solid line) from Eq. (|4.17 ) and as obtained by keeping only the first two terms of the sum 



(dashed line). The typical susceptibility from Eq. ( 4.18|) is represented by the dotted horizontal 
line. 



FIG. 4. First harmonic of the persistent current in a ring with r = h/ a = 0.9 as a function 
of kpa (solid line) for a cut-off length = 6a according to Eqs. ( 4.27 )-( 4.28| ) together with the 
contribution coming from type-I trajectories (dashed line). The typical persistent current from 
Eq. ( [4.29| ) is represented by the dotted horizontal line. 



FIG. 5. First harmonic of the typical persistent current in rings of different thickness (r = h/a) 
for various cut-off lengths Lc = 30a (circles), 6a (diamonds) and 3a (triangles) according to 
Eq. ( [4.29| ). Filled symbols correspond to the total persistent current and lay approximately on 
a horizontal line for each Lc, consistent with the asymptotic behavior of Eq. ( 4.391 ) indicated 
by arrows on the extreme right of the plot. Unfilled symbols represent the contributions from 
each type of trajectory and are joined by dotted lines (type-I) and dashed lines (type-II). These 
guide-to-the-eye exhibits the approximate behavior of Eqs. ( |4.3^ ) and shows how the r character- 
istic of the switching from one type of trajectories to the other increases with temperature. 



FIG. 6. a) Trajectory from the family (1,1) of the square billiard. The abscissa xq of the 
intersection of the trajectory with the lower side of the square, together with the label e = ±1 
precising the sense of the motion, label the trajectories inside the family, b) Trajectory from the 
family (2,1) of the square billiard, illustrating the fiux cancelation occurring for other periodic 
trajectories than those in the (1,1) family (or their repetitions). 
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FIG. 7. The method of images: The Green function G{q, q') is constructed from the free Green 
function by placing a source point at each mirror image Qi of the actual source q. To each of 
the resulting contribution G^{qi,q') is associated a classical trajectory (solid line). This latter is 
obtained from the straight line joining qi to q' (dash line) by mapping all its intersected images 
back onto the original billiard. 



FIG. 8. (a) Magnetic susceptibility of a square as a function of kpa from numerical calculations 
(dotted line) at zero field and at a temperature equal to 10 level-spacings. The number of electrons 
is = gs{kpa)'^ / {^tt). The full line shows our semiclassical approximation (Eq. (|5.1^ )) taking into 
account only the family (11) of shortest orbits with the temperature correction factor Rt{Lii). The 
period vr/\/2 of the quantum result indicates the dominance of the shortest periodic orbits enclosing 
non-zero area with length Ln = 2^/2a. (b) Susceptibility x as a function of the normalized flux 
through the sample (at a Fermi energy corresponding to ~ 400 enclosed electrons) from Eq. ( ^.18 ) 
(solid) and numerics (dotted). The susceptibility arising from the stationary-phase integration 
(Eq. ( |5.13D ) shown as the dashed line diverges at 99 — > 0. 



FIG. 9. (a) Average magnetic susceptibility of an ensemble of squares differing in size as 
a function of kpa for various temperatures (8, 6 and 4 level spacings for the three triplets of 
curves from below) and a flux (p = 0.15. Solid line: average of the semiclassical approximation 
to x*-^^ according to the analytical result of Eq. ( |5.21| ). Dotted line: average of x^^^ obtained by 
using Eq. (2.20c) and exact diagonalization. Dashed curve: average of the canonical susceptibility 
calculated directly from Eq. ( ^.33 ) after the exact diagonalization. The considerable agreement 
between the solid and dotted curves illustrates the precision of the semiclassical apprroximation, 
while the agreement between the dotted and dashed lines shows the applicability of the thermo- 
dynamical expansion Eq. ( 2.20[) . (b) Flux dependence of the averaged susceptibility normalized to 
Xn = XikpaRj^^Lii) at kpa ~ 70 from the semiclassical expression Eq. (|]2lD (sohd) and numer- 
ical calculations (dashed). The thick solid (dashed) curve denotes an average of the semiclassical 
(numerical) result over an ensemble with a large dispersion of sizes which is denoted by (x) (see 
text). The shift of the numerical with respect to the semiclassical results reflects the Landau sus- 
ceptibility (due to in Eq. ( 2.19| )) and effects from bouncing-ball orbits (see section VII A| ) not 
included in the semiclassical trace. 



FIG. 10. Solid: low temperature limit of the average susceptibility x*^^-* of an ensemble of 



squares, as given by Eq. (|5.27D and normalized by xa^ as defined in the previous figure caption. 
Dashed: contribution of the family (1,1) to this result. Even in the very low temperature regime 
the magnetic response is dominated by the (1,1) family except for the singularity which develops 
at zero magnetic field. 



FIG. 11. Solid: exact (Presnel) function C{ip) as given by Eq. ( 5.11| ). Dashed: approximation 
of C{^p) by the Bessel function Jo{32ip/tt'^) (see text). 
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FIG. 12. Solid: function F(^) (see Eq. (|6.29|) ) describing the magnetic field dependence of the 
average susceptibility for an ensemble of chaotic microstructures. Dashed: quadratic approximation 
(^V6-(3C(3)+2^Vl5)e')ofF(e). 



FIG. 13. Grand-canonical susceptibility of a square potential well as a function of magnetic 
flux Lp = Ha^/^Q. The full lines always denote the quantum mechanical results. Panel (a): x/xl 
calculated at a Fermi energy of 2140 enclosed electrons at a temperature kT/A = 8. Dashed 
(dotted) line: Semiclassical result due to bouncing-ball orbits from Eq. ( 7.11| ) with action ac- 
cording to the exact expression of Eq. ( |7.3|) , (quadratic approximation Eq. (|7.12| )). Dashed-dotted 
line: Susceptibility contribution from family (11) from Eq. (|5.18D with offset of —80 for reasons of 
representation, (b) Dashed line: Semiclassical contribution (Eq. ( |7.13| )) from bouncing-ball orbits 
for 1440 electrons and kT/A = 7. The lower value of kp makes it necessary to describe the actions 
by Eq. ([7.3|). (c) same as in (b) but for a low temperature kT/A = 2 for which repetitions are 
important and the use of Eq. (f7l3|) is necessary to approach the quantum results. 



FIG. 14. Schematic representation of a typical flux-enclosing bouncing-ball orbit with cy- 
clotron radius re- The dashed lines denote the limits of the i?-dependent range of bouncing-ball 
orbits. 



FIG. 15. De Haas-van Alphen like oscillations of the susceptibility of a square at magnetic 
fluxes corresponding to Tc < a/2 for 2140 electrons at kT/A = 8. Full line: quantum calculations; 
dashed line: analytical semiclassical result from cyclotron orbits according to Eq. ( [7.19| ). 

FIG. 16. Contour of integration in the complex energy plane used to evaluate the integral 
Eq. (|Al| ). The derivative f'{E — n) of the Fermi function has a periodicity of 2i7r//? and double 
poles located at En = fJ- + i{2n + l)vr//3 (n being a positive or negative integer). Moreover, at 
low temperature, f'{E — //) is essentially zero except for a narrow band of width (3~^ near the line 
Re(i?) = //. With this contour of integration, the integrand of Eq. ( |A1| ) has to be evaluated only 
in the small domain [fj, — /3~^, /x -|- x [0, 2i7r/3~^] where a linearized approximation of the action 
is accurate. 



FIG. 17. Integration paths on the invariant tori used to compute the action integrals Ii and I2 
for the circular and the ring geometries, a) Path Ci (thick dashed) and C2 (thick solid) for the circle 
and type-I trajectories of the ring. The straight part of C2 is along a classical trajectory (thin solid), 
whose caustic (dotted) is outside the inner disk, b) Path C2 (thick solid) for type-II trajectories of 
the ring. Path Ci is similar as in a) and therefore not shown. The straight parts of C2 are along 
classical trajectories. We indicated one of them (solid thin) and its caustic (dotted) laying inside 
the inner disk. The thick-dashed line joining the straight parts of C2 is a guide-to-the-eye putting 
in evidence the simple form of Eq. (p3|). 
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TABLES 





D°^''{E)/D 


x/xl 


x/xl 


chaotic 






{k^af 


regular 




{k^afl^ 


{kpo) 



TABLE I. {kpa) dependence of the oscillating part of the density of states and of the magnetic 
response depending on the absence (chaotic case) or the presence (regular case) of continuous 
families of periodic orbits for two-dimensional billiard-like microstructures. {D = (gsinA) / (Q.irh'^) 
is independent of the nature of the dynamics.) 
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